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ABSTRACT 

The late time optical and near-IR line profiles of many core-collapse supernovae 
exhibit a red-blue asymmetry as a result of greater extinction by internal dust of 
radiation emitted from the receding parts of the supernova ejecta. We present here 
a new code, DAMOCLES, that models the effects of dust on the line profiles of 
core-collapse supernovae in order to determine newly formed dust masses. We find 
that late-time dust-affected line profiles may exhibit an extended red scattering wing 
(as noted by Lucy et al. (1989)) and that they need not be flux-biased towards the 
blue, although the profile peak will always be blueshifted. We have collated optical 
spectra of SN 1987A from a variety of archival sources and have modelled the Ha 
line from days 714 to 3604 and the [O i] 6300,6363 A doublet between days 714 and 
1478. Our line profile fits rule out day 714 dust masses > 3 x 10 -3 M 0 for all grain 
types apart from pure magnesium silicates, for which no more than 0.07 M 0 can be 
accommodated. Large grain radii (^ 0.6 pm) are generally required to fit the line 
profiles even at the earlier epochs. We find that a large dust mass (> 0.1 M 0 ) had 
formed by day 3604 and infer that the majority of the present dust mass must have 
formed after this epoch. Our findings agree with recent estimates from spectral energy 
distribution fits for the dust mass evolution of SN 1987A and support the inference 
that the majority of SN 1987A’s dust formed many years after the initial explosion. 

Key words: radiative transfer - supernovae: general - supernovae: individual: SN 
1987A - ISM: supernova remnants. 


1 INTRODUCTION 

Core-collapse supernovae (CCSNe) have long been thought 
to be potential dust factories (Hoyle & Wickramasinghe 
1970; Kozasa et al. 1991; Todini & Ferrara 2001). However 
over the previous decade observations at mid-infrared (mid- 
IR) wavelengths of warm dust emission from CCSNe had 
suggested that the quantities of dust produced, typically < 
10 -3 M 0 during the first 1000 d (Sugerman et al. 2006; 
Meikle et al. 2007; Kotak et al. 2009; Andrews et al. 2010; 
Fabbri et al. 2011) were much less than the 0.1-1.0 M 0 of 
dust per CCSN estimated to be needed (Morgan & Edmunds 
2003; Dwek et al. 2007) in order to account for the very 
large dust masses measured in some high redshift galaxies 
(Omont et al. 2001; Bertoldi et al. 2003; Watson et al. 2015). 
However, recent Herschel far-IR and sub-mm observations 
of cold dust masses as high as 0.2-0.8M 0 in several young 
supernova remnants have resulted in a re-evaluation of the 
rate of dust production by CCSNe (Barlow et al. 2010; Mat- 
suura et al. 2011; Gomez et al. 2012). The Hershel dust mass 
estimates were based on fitting dust spectral energy distribu¬ 


tions (SEDs) that peaked at far-IR wavelengths. Following 
the end of the Herschel mission in 2013 there is likely to 
be a long wait for far-IR facilities with comparable or better 
sensitivities than Herschel to become available, providing an 
incentive to make use of alternative methods to estimate the 
dust masses that form in supernova ejecta. 

The absorption and scattering of optical or near-IR ra¬ 
diation by newly-formed dust within the ejecta of super¬ 
novae can result in an asymmetry between the red and 
blueshifted components, with redwards emission from the 
far side of the ejecta undergoing greater absorption. Lucy 
et al. (1989) identified a progressive blueshifting of the 
[O i] AA6300,6363 A doublet from SN 1987A between days 
529 and 739 after outburst, with the doublet in the later 
spectrum being blueshifted by ~ 600 km s -1 . Since then, 
such red-blue asymmetries have been frequently observed 
in the late-time (> 400 d) spectra of supernova ejecta and 
there is now a growing data base of such observations (e.g. 
Lucy et al. (1989); Fabbri et al. (2011); Mauerhan & Smith 
(2012); Milisavljevic et al. (2012)). 

SN 1987A as an archetypal object is critical to our grow- 
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Figure 1. Archival data showing the evolution of the Ha and 
[O i] line profiles from SN 1987A at the earlier of the epochs 
considered. The spectral gaps at the last two epochs correspond 
to where narrow line emission from the equatorial ring has been 
removed. The spectra have been continuum-subtracted and offsets 
have ben applied for display purposes. 


ing understanding of the formation and evolution of dust in 
CCSNe. Since its outburst there have been numerous obser¬ 
vations at many wavelengths and many epochs. Mid-infrared 
emission from warm dust (T~400 K) was observed by day 
~450 (Roche et al. 1989; Bouchet et al. 1991; Wooden et al. 
1993) and by day 775 the emitting dust mass was estimated 
to have been between ~ 5 — 20 x 10 -4 M@ (Wooden et al. 
1993; Ercolano et al. 2007; Wesson et al. 2015). Beginning 
from 23 yr after outburst, the Herschel Space Observatory 
detected much larger quantities (0.4-0.8 Mq of T~20 K 
cold dust emitting at far-IR and submillimetre wavelengths 
(Matsuura et al. 2011, 2015). This emission has been con¬ 
firmed by ALMA observations to originate from the ejecta 
of SN 1987A (Indebetouw et al. 2014). 

We here seek to model the effects of dust on line pro¬ 
files with a view to providing both an alternative way of de¬ 
termining dust masses formed in the ejecta of CCSNe and 
in order to investigate the effects of dust on the shapes of 
line profiles emitted from these objects. We present a new 
code, DAMOCLES (Dust Affected Models Of Character¬ 
istic Line Emission in Supernovae), that utilizes a Monte 
Carlo methodology in order to model line profiles in ex¬ 
panding atmospheres. The code can treat dust composed 
of multiple species and grain sizes with variable ejecta den¬ 
sity and velocity distributions. Both clumped and smooth 
geometries may be modelled. 

In this paper, we collate optical spectra from the 
archives of four different telescopes in order to study the 
effects of dust formation on the Ha line and on the 
[O i] A6300,6363 A doublet. We model epochs spanning a 
range of approximately 8 yrs from the first indications of 
blueshifting in the Ha line between days 600-700, using both 
smooth and clumped geometries. We compare our derived 
dust masses to those obtained by Wesson et al. (2015, here¬ 
after 2015) and Dwek & Arendt (2015, hereafter DA15) and 


Figure 2. Archival data showing the evolution of the Ha line 
profile from SN 1987A at the later epochs. The spectral gaps 
correspond to where narrow line emission from the equatorial ring 
has been removed. The spectra have been continuum-subtracted 
and offsets applied for display purposes. 

consider the implied dust formation rate. We present our 
testing of the new code against analytical cases and previ¬ 
ously published optically thick models (Lucy et al. 1989). 
We also investigate the sensitivity of line profiles to each of 
the variables and note the range of signatures that observed 
line profiles may exhibit in the presence of dust. 

In Section 2, we detail the observed spectra that we 
used for our modelling. In Section 3, we discuss the details 
of the DAMOCLES code and in Section 4 we present our 
testing of the code and our parameter sensitivity analyses. 
Our modelling of the Ha and [O i] A6300,6363 A lines is 
presented in Section 5 and we discuss our findings in Section 
6 . 


2 ARCHIVAL SPECTRA OF SN 1987A 

SN 1987A has been the most intensively observed supernova 
in history, with a wealth of both spectral and photometric 
data available to model. From the archives of a number of 
different telescopes we have collated optical spectra acquired 
over a wide range of epochs. At the earlier epochs we use 
spectra obtained by the Anglo-Australian Telescope (AAT) 
and the Cerro Tololo Inter-American Observatory (CTIO) 
and at later epochs we use spectra from the archives of the 
Hubble Space Telescope ( HST) and the Very Large Telescope 
(VLT). An explosion date of 1987 February 23 is adopted 
throughout and epochs are measured relative to this date. 
Full details of all observations may be found in Table 1. The 
spectral resolutions of the grating spectrograph observations 
are listed in column 7, while column 8 lists the spectral 
resolving powers of the echelle spectrograph observations. 

Wavelength ranges encompassing the Ha line and 
[O i] AA6300,6363 A doublet were selected in order to trace 
their evolution from day 524, near the time of the first indi¬ 
cations of dust formation (Wooden et al. 1993), to day 8020, 
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Table 1. Details of the archival data for SN 1987A. 


Date 

Age 

(d) 

Telescope 

Inst 

^min 

(A) 

^max 

(A) 

Res. 

(A) 

Res. Power 

Reference 

31 Jul 1988 

524 

AAT 

FORS 

5500 

10190 

20 


Spyromilio et al. (1991) 

26 Oct 1988 

611 

AAT 

UCLES 

6011 

7336 


30000 

Hanuschik et al. (1993); Spyromilio et al. (1993) 

27 Dec 1988 

673 

AAT 

UCLES 

5702 

10190 


30000 

Hanuschik et al. (1993); Spyromilio et al. (1993) 

06 Feb 1989 

714 

CTIO-1.5m 

Cass. 

6420 

10380 

16 


Phillips et al. (1990) 

09 May 1989 

806 

CTIO-1.5m 

Cass. 

6430 

10330 

16 


Phillips et al. (1990) 

12 Jan 1990 

1054 

CTIO-4m 

RC 

3565 

10000 

11 


Suntzeff et al. (1991) 

12 Mar 1991 

1478 

CTIO-4m 

RC 

3245 

9175 

11 



30 Mar 1992 

1862 

HST 

STIS 

4569 

6818 

4.4 


Wang et al. (1996) 

14 Mar 1993 

2211 

HST 

STIS 

4569 

6818 

4.4 


Wang et al. (1996) 

07 Jan 1995 

2875 

HST 

STIS 

4569 

6818 

4.4 


Chugai et al. (1997) 

23 Sep 1996 

3500 

HST 

STIS 

4569 

6818 

4.4 



05 Jan 1997 

3604 

HST 

STIS 

4569 

6818 

4.4 



10 Dec 2000 

5039 

VLT 

UVES 

4760 

6840 


50000 

Groningsson et al. (2006, 2007) 

06 Oct 2002 

5704 

VLT 

UVES 

4760 

6840 


50000 

Groningsson et al. (2006, 2007, 2008) 

21 Mar 2005 

6601 

VLT 

UVES 

4760 

6840 


50000 

Groningsson et al. (2006, 2007) 

23 Oct 2007 

7547 

VLT 

UVES 

4760 

6840 


50000 

Groningsson et al. (2007) 

07 Feb 2009 

8020 

VLT 

UVES 

4800 

6800 


50000 

Tziamtzis et al. (2010) 


near the current era. Optical spectroscopy obtained at the 
AAT using the Faint Object Red Spectrogaph (FORS) dur¬ 
ing the first 2 yr after outburst was kindly supplied by Dr 
Raylee Stathakis (Spyromilio et al. 1991; Hanuschik et al. 
1993; Spyromilio et al. 1993) and optical spectra from the 
CTIO were donated by Dr Mark Phillips (Suntzeff et al. 
1991). 

The evolution of the Ha and [O i] line profiles is pre¬ 
sented in Figs 1 and 2. At later epochs, the broad Ha profile 
emitted by the ejecta becomes contaminated by narrow line 
emission from the equatorial ring. These lines have been re¬ 
moved for the purposes of modelling the broad line. A con¬ 
tinuum fit has been subtracted from each spectrum and a 
velocity correction has been applied for a recession velocity 
of 287 km s -1 (Groningsson et al. 2008). 


2.1 Contamination of the Ha profiles 

The Ha profile at day 714 exhibits a very slight inflection 
visible at V « +900 km s -1 . By day 806, this slight inflection 
has developed into a noticeable shoulder in the line profile 
of Ha (see Fig. 11). 

Although these features are similar in nature to features 
produced by dust absorption in the flat-topped region (as 
discussed in Section 4.3.5), we conclude that this shoulder 
is an early appearance of the unresolved [N n] A6583 A line 
from the equatorial ring (Kozma & Fransson 1998b). Unre¬ 
solved nebular [N ii] lines at A = 6583 A and A = 6548 A 
either side of the Ha rest-frame velocity at 6563 A are cer¬ 
tainly seen by day 1054 and have to be removed in order 
to consider the evolution of the broad Ha profile (see Fig. 
1). We do not remove this potential contaminant at earlier 
epochs but try to fit the broad line profiles around it. 

By day 1054, all three of the narrow nebular lines are 
strong. They remain unresolved in the low spectral resolu¬ 
tion CTIO data at days 1054 and 1478 and therefore contam¬ 
inate the entire central region of the Ha line profile. Their 
presence renders two CTIO Ha profiles from days 1054 and 
1478 unusable for modelling purposes. The HST and VLT 


Table 2. Ha FWHM and the HWZI determined by the zero in¬ 
tensity velocity on the blue side of the line. The tabulated line 
widths have been corrected for the relevant instrumental resolu¬ 
tion. 


day FWHM (A) HWZI (A) 


524 

3200 

3600 

611 

2700 

3400 

673 

1600 

3700 

714 

3100 

4500 

806 

3200 

5500 

1054 

2100 

5600 

1478 

1400 

6600 

1862 

1600 

6800 

2211 

1400 

6700 

2875 

2700 

6700 

3500 

3500 

7000 

3604 

2100 

7000 


Ha profiles at later epochs (> 1862 d) have a higher spectral 
resolution and it was therefore easier to remove the narrower 
[N n] and Ha lines from the broad Ha profiles (for example 
Figs 1 and 2). Although this does remove a potentially in¬ 
formative section of the profile (+500 km s^ 1 < v < +1500 
km s -1 ), we achieve good fits to the overall line profiles at 
these epochs. 

2.2 The evolution of the maximum and minimum 
velocities 

For a freely expanding medium, the velocity of any fractional 
radial element should not change with time. The maximum 
velocity of any line-emitting region is therefore expected to 
be constant. However, at the epochs we consider here, it 
appears that the maximum velocities of the Ha line, as de¬ 
termined by the velocity at zero intensity on the blue side, 
generally increase over time (see Table 2). We attribute this 
to the start of the freeze-out phase in the outer regions of the 
ejecta, while the hydrogen neutral fraction is still increasing 
in the denser inner regions (Danziger et al. 1991; Fransson 
& Kozma 1993). 
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The onset of a fixed ionization structure in the ejecta 
causes the rate of Ha flux decline to slow. Since the outer, 
faster moving regions reach this state at earlier times than 
the inner, slower moving regions, the relative flux contri¬ 
bution of the outer regions is increased. At early epochs 
(t < 900 d) the flux contribution from hydrogen in the 
core dominates the overall Ha flux, whereas at later epochs 
(t. > 900 d) the flux from the envelope dominates (Fransson 
& Kozma 1993; Kozma & Fransson 1998a). This shift likely 
explains apparent broadening of the line with the higher ve¬ 
locity material becoming increasingly noticeable in the line 
profiles. This may also explain the increase in half-width 
zero intensity (HWZI) velocities at these epochs with the 
relative flux from the very densest regions dropping more 
rapidly relative to the outer line-emitting region. The full 
width at half-maximum (FWHM) remains relatively steady 
(see Table 2). However, the FWHM values presented in Ta¬ 
ble 2 were difficult to determine accurately since the peak of 
the broad line profile is contaminated by narrow line emis¬ 
sion from the equatorial ring. 


3 THE DAMOCLES CODE 

Monte Carlo methods have long been used to model radia¬ 
tive transfer problems in diverse environments and there are 
several examples of codes which apply the technique in appli¬ 
cation to supernovae (for example Maeda et al. (2003); Lucy 
(2005); Jerkstrand et al. (2012); Owen & Barlow (2015)). 
Whilst there are numerous codes that treat dust or gas or 
both in order to produce an overall SED, there is a dearth 
of codes designed to focus on the shapes of individual line 
profiles. Although a velocity field is naturally considered in 
codes that seek to reproduce the spectra of supernovae, ab¬ 
sorption and scattering by dust is not and thus the result¬ 
ing shapes of line profiles are potentially unrepresentative of 
those emerging from dusty ejecta at late times. 

In this work we aim to model single or doublet line pro¬ 
files produced by a moving atmosphere in a dusty medium. 
Since a comparatively small wavelength range is consid¬ 
ered, a fully self-consistent radiative transfer model is un¬ 
necessarily expensive. Instead any energy packet that is ab¬ 
sorbed during the simulation may simply be removed on the 
grounds that it would be reemitted outside the wavelength 
range of interest. This approach is clearly not applicable to 
SED radiative transfer models that treat continuum emis¬ 
sion from dust. The extinction due to dust is assumed to 
be temperature-independent and it is therefore unnecessary 
to iteratively calculate the temperature of the ejecta as in a 
fully self-consistent calculation of the SED. Though clearly 
the total energy transferred through the medium is not con¬ 
served in the wavelength range of interest, the signature of 
the normalized line profile is preserved. 

The DAMOCLES code builds on the work of Lucy 
et al. (1989) who employed a similar approach to model 
the broad [O i] A6300,6363 A doublet seen in SN 1987A at 
early epochs (up to ~ day 775). It models the transport of 
initially monochromatic energy packets through a smooth 
or clumped dusty medium having a smooth velocity field. 
The velocity field and the inner and outer ejecta radii are 
free parameters. The late-time (> 400 d) line emission is 
assumed to be optically thin, with an emissivity distribu¬ 


tion proportional to the square of the local gas density, i.e. 
proportional to the product of the recombining proton and 
electron densities in the case of Ha or to the product of the 
neutral oxygen and electron densities in the case of collision- 
ally excited [O i] emission. 

3.1 The energy packet formalism 

The initial radiation field is inherently tied to the distribu¬ 
tion of gas throughout the supernova ejecta which is declared 
as a power law p(r) oc r -/3 between Ri n and Rout- Rout is 
calculated directly from the epoch of the line to be modelled 
and the declared maximum line velocity. The emissivity dis¬ 
tribution is also specified as a power law with i(p) oc p k . 
However this is generally taken to be i(r) oc r~ 213 since the 
majority of lines modelled are optically thin recombination 
lines or collisionally excited lines and therefore i(p) oc p 2 . 
The radiation is quantised into monochromatic packets with 
equal energy Eq = nhi/o- In Monte Carlo simulations (that 
model non-moving media) packets are usually taken to be of 
constant energy. When the frequency of a packet is altered 
after an event, the energy of that packet is kept constant and 
the number of real photons contained within it is assumed 
to change. However, in the case of dust scattering, the num¬ 
ber of real photons is conserved and thus the energy of the 
packet is altered. This is most easily achieved by weighting 
each packet over all scattering events as 


scat 

where w is the weight of the packet and v and v' are the 
frequencies of the packet before and after the scattering 
event respectively. The final energy of each packet is then 
E = wEo, where Eo is the initial energy of the packet. 

The emissivity distribution is calculated by dividing the 
ejecta into a specifiable number of shells between 77„ and 
Rout overlaid on the Cartesian grid and the number of pack¬ 
ets to be emitted isotropically in each shell calculated ac¬ 
cording to the specified emissivity and density power laws. 
For each packet a location within that shell and an initial 
trajectory is randomly sampled from an isotropic distribu¬ 
tion such that 

cj> = 27 rq ( 2 ) 

cos 0 = 2£ — 1 (3) 

where 0 < r] < 1 and 0 < £ < 1 are random numbers, <j> is 
the azimuthal angle and cos 9 is the radial direction cosine. 
At emission and at each scattering event the frequency of 
the packet is recalculated according to the specified radial 
velocity field v(r) oc V r „a. x r a (see Section 3.3). 

3.2 The geometry of the ejecta and the grid 

The supernova ejecta is approximated by a three- 
dimensional Cartesian grid, each cell of which is assumed to 
have uniform density and composition. The grid is a cube 
with sides of width 2R out and a declarable number of divi¬ 
sions. After the initial emission of energy packets, the gas 
plays no further role in the simulation and thus only dust 
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Figure 3. Red: benchmark models for optically thin (r = 0) line profiles with fractional velocity v ex r. Left to right: initial emissivity 
profiles i(r ) oc r~ with f} = 0.0,1.0 and 2.0. Cases with Ri n /R 0 ut = 0.2 are on the top and with Ri n /R 0 ut = 0.0 on the bottom. The 
presence of a plateau in the upper plots is due to the finite inner radius (detached shell). Blue: the analytical case with i(u ) ~ 1 — u 211 
except in the case of j3 = 1 where i(u) ~ — log u. Peak fluxes are scaled to unity. 


properties are considered. By default, the dust is coupled 
to the gas (although it may be decoupled) and thus follows 
the smooth distribution described above (p oc r ~ 18 ). The 
dust density in each cell is therefore calculated accordingly 
and any cell whose centre falls outside of the bounds of the 
supernova ejecta has density set to zero. 

It is worth noting that if a constant mass-loss rate is re¬ 
quired, the exponent of the velocity profile and the exponent 
of the density profile are not independent. A constant mass 
loss rate implies that 4npvr 2 oc k , where k is a constant, and 
thus for v oc r a and p oc r - ' 3 , we require that /3 — a = 2. 
However, it is possible that the supernova event may have 
induced a mass-flow rate that is not constant with radius 
and thus both exponents may be declared independently. 

It is known from SED modelling that clumped envi¬ 
ronments produce very different results to environments as¬ 
sumed to have a smooth distribution of dust and gas. Specif¬ 
ically, clumped models tend to require a higher dust mass 
in order to reproduce a similar level of infrared dust emis¬ 
sion compared to a smoothly distributed model. The ca¬ 
pacity for modelling a clumped dusty medium is therefore 
included in the code. The fraction of the dust mass that is in 
clumps is declared (mf rac ) and the total volume filling factor 
of the clumps (/) is also specified. Dust that is not located 
in clumps is distributed according to a smooth radial pro¬ 
file. The clumps occupy a single grid cell and their size can 
therefore be varied by altering the number of divisions in 
the grid. The clumps are distributed stochastically with the 
probability of a given cell being a clump proportional to the 


smooth density profile (i.e. p(r) oc r ^). The density within 
all clumps is constant and is calculated as 

_ Afdumpa mfrac A/tot / . \ 

Pclump “ w - iMRLt - Rfn) 

where Mtot is the total dust mass, M c i umpa is the total dust 
mass in clumps and Uciumps is the total volume occupied by 
clumps, mfrac and / are defined as above. 

3.3 The radiative transport mechanism 

Following emission, a packet must be propagated through 
the grid until it escapes the outer bound of the ejecta at 
.Rout. The probability that the packet travels a distance l 
without interacting is p(l) = e~ nal = e -T where n is the 
grain number density, a is the cross-section for interaction 
and r = nal for constant n and a (as in a grid cell). Not¬ 
ing that the probability that a packet will interact within a 
distance Z is 1 — e -T , we may sample from the cumulative 
probability distribution to give: 

i = 1 - e -T => r = -ln(l - £) (5) 

where 0 < £ < 1 is a random number sampled from a uni¬ 
form distribution. The frequency of the photon packet and 
the mass density of the cell are then used to calculate the 
opacity of that cell. Using the fact that na = ftp, the dis¬ 
tance l that the packet travels before its next interaction is 
calculated. If this value is greater than the distance from 
its position to the edge of the cell then the packet is moved 
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along its current trajectory to the cell boundary and the 
process is repeated. If the distance is less than the distance 
to the boundary then an event occurs and the packet is ei¬ 
ther scattered or absorbed, with the probability of scattering 
equal to the albedo of the cell 


If the packet is absorbed then it is simply removed from 
the simulation as discussed above. If the packet is scattered 
then a new trajectory is sampled from an isotropic distri¬ 
bution in the comoving frame of the dust grain and the 
frequency of the packet is recalculated using Lorentz trans¬ 
forms subject to the velocity at the radius of the interaction 
(see Appendix A for further details). This process is repeated 
until the packet has either escaped the outer boundary of 
the supernova ejecta or has been absorbed. Escaped pho¬ 
ton packets are added to frequency bins, weighted by w, in 
order to produce an overall emergent line profile. For our 
grain radii between 0.35/rm and 3.5/rm, Mie theory predicts 
grains to be forward scattering (g = 0.7 — 0.8). However, 
models that incorporate these values with the Henyey & 
Greenstein (1941) phase function are not significantly dif¬ 
ferent to isotropically scattering grain models. 


as described above. Including the electron scattering mech¬ 
anism in models is optional. 

In the majority of cases the electron scattering optical 
depths are not high enough to discernibly affect the overall 
shape of the profile. However, there may be some early epoch 
cases (the concept is discussed for SN 2010jl by Fransson 
et al. (2014)) where the electron scattering optical depths 
are high enough to have a significant effect on the observed 
profiles. 


3.5 Technical details 

DAMOCLES is written in FORTRAN 95 and parallelized 
for shared memory machines using OPENMP. It has been 
developed on and currently runs on a MacBook Pro 11.2 
quad core with Intel Core i7 2.8GHz processors and 16GB 
of memory. A typical, medium resolution simulation using 
125,000 grid cells and 10 5 packets takes approximately 15 
s to run. The number of packets transported and the total 
dust optical depth are the most important factors in deter¬ 
mining runtime. We intend to make DAMOCLES available 
in the public domain after some further development and 
documentation. 


3.4 Properties of the Dusty Medium 

Dust of any composition for which optical data are avail¬ 
able may be used and the relative abundances of the species 
may be declared by the user. A grain size may be specified 
for each species. The extinction due to dust is only depen¬ 
dent on the cross-sectional area of the grains and not on the 
overall distribution. It is therefore not possible to determine 
details of a grain size distribution and only to constrain a 
single grain size parameter. The capacity to declare a size 
distribution is however included for the sake of ease of com¬ 
parison with SED models. Mie theory and optical proper¬ 
ties are used to calculate the overall Qnhs{ v ) an d QscaM for 
each species and the derived opacities are summed over each 
species weighted according to their relative abundances. 

As will be discussed in Section 4, the effects of scatter¬ 
ing on the shapes of line profiles can potentially be quite 
pronounced and it is therefore important to consider the 
effects of electron scattering as well as those of dust scat¬ 
tering. Electron densities are estimated from the observed 
luminosity of Ho- using an average temperature of 10,000 K 
and assuming that the electron density distribution is cou¬ 
pled to the emissivity distribution. The total optical depth 
to electron scattering between Ri n and R ou t is then calcu¬ 
lated. Electron scattering is treated in an identical manner 
to dust scattering, with r = ra U st + r e in each cell. If, for a 
given packet, an event occurs, it is first calculated whether 
this is an electron scattering event or a dust event (either 
scattering or absorption) by considering the ratio of the op¬ 
tical depths to each. If the packet is scattered by an electron 
then the velocity of that electron is calculated by considering 
the bulk velocity at that radius and adding a thermal veloc¬ 
ity component following the formalism described by Hillier 
(1991). The scattering process is then identical to that for 
dust. If the event is a dust event then the process continues 


4 COMPARISON OF DAMOCLES MODELS 
WITH ANALYTICAL AND PREVIOUSLY 
PUBLISHED RESULTS 

There is a general lack of published models in the literature 
that consider dust absorption-affected asymmetric line pro¬ 
files. We therefore tested the code by comparing the results 
to optically thin profiles that may be derived analytically. 
We then tested the absorption and scattering components 
of the code by comparing our results for the case of an opti¬ 
cally thick medium with those derived by Lucy et al. (1989) 
in their Model II and Model III scenarios. 


4.1 Comparison of DAMOCLES models with 
analytical results 

Analytical profiles may be calculated in the dust-free case. 
We ran a number of models based on the methods of Gerasi- 
movic (1933) who derived equations for line profiles emitted 
from a transparent expanding shell. 

Describing the fractional expansion velocity of the shell 
as v(r) oc r“ with a^0 such that v(r) = y r where V(r) 
and Umax represent physical velocities and u max = 1, the 
energy emitted by the nebula between line of sight velocities 
u and u + du is proportional to 

J i(r)r sin($) r dd dr (7) 

where i(r) represents the emission per unit volume at radius 
r and 9 is the angle to the observer’s line of sight. We adopt 
inner radius Ri n = q and outer radius R ou t = 1 such that 
q — Rin / Rout- 

Setting i(r) oc r~ 2/3 (for a recombination or collisionally 
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Figure 4. Benchmark models for line profiles with v oc r, i(r) oc constant and a filled sphere with Ri n /R 0 ut = 0. Pure dust absorption 
models (ui = 0) are presented in the left-hand plots, whilst partially scattering models are presented on the right ( uj = 0.6) as per Lucy 
et al. (1989) Models II and III. All resulting profiles have been scaled to unit flux at their peaks. 


excited line emitted from a medium with an assumed density 
profile for the emitter p oc r~ 13 ) then gives 


i(u) du ' 


du 

2/3-3+a 

IU i 

du 

28-3 + a 


rS i 

Je o 


2,8-3 

cos “ 9 sin 9 d 9 


9 


2/3 - 3 + a 


( 8 ) 


for 2 ^~ 3 —1 where i(u) du is the energy emitted in a 

volume element and 9o and 9\ are the bounds of this element. 
The case 2|3 q ) 3 = — 1 results in a logarithmic relationship. 

In the case of a “filled” nebula, i.e. one where the inner 
radius is vanishingly small in comparison to the outer radius, 
we obtain 


i{u) du ' 


du 


(2/3 — 3 + a)u~ 


mi 


1 — u 


2/3 —3 + o 


) (9) 


If the nebula is not “filled”, that is to say, the inner 
radius is some fraction of the outer radius and the remnant 
is a detached shell, the above formula becomes valid only 
from some critical value u! = q a to u = 1. For u < u', we 
obtain 


i ( u ) du ~ ± { 2/3 _3 + Q) (^~ 1 ) (10) 

and therefore the top of the line is flat while the sides are 
sloping. 

Crucially, the width of the flat section is determined by 
u! = q a or simply u! = q in the case where u oc r, whilst the 
shape of the profile outside of the flat top is described by 
equation 9. 

Profiles with a variety of shapes may be derived from 
these formulae depending on the relative values of a and /3. 
Here we consider three main families of curves: 


(ii) i(u) ~ 1 — u 7 (a > 0, 2/3 — 3 + a < 0) 

(iii) i( u) ~ — log u (a > 0, 2/3 — 3 + a = 0) 

where 7 is defined as 7 = | 2l3 ~^+ a |. 

Models are presented for each of these cases, both for a 
filled nebula and for a shell structure with Ri n /R 0 ut = 0.2. A 
velocity profile v oc r appropriate for supernova ejecta in the 
free expansion phase is used throughout (Li & McCray 1992; 
Xu et al. 1992; McCray 1996; Baron et al. 2005). Values of 
/3 = 0,1 and 2 are adopted. Fig. 3 illustrates the excellent 
agreement between the analytical case and the models. All 
fluxes are scaled to unity at the peak. 


4.2 Comparison of DAMOCLES models with 
previously published results 

In addition to the tests for optically thin lines described 
above, we also compared our outputs to those derived by 
Lucy et al. (1989) in order to assess the accuracy of the 
scattering and absorption aspects of the code. We consider 
two similar cases, equivalent to Models II and III of Lucy 
et al. (1989). In the first case, dust with zero albedo (pure 
absorption) is uniformly distributed throughout a filled neb¬ 
ula with a velocity profile v oc r. In the second case, the same 
scenario is considered but in a medium of dust with albedo 
uj = 0 . 6 . 

In the first case, the profile may once again be derived 
analytically from the basic geometry using the fact that ra¬ 
diation will be attenuated by a factor e -2 ” 1 '” between points 
with line-of-sight fractional velocities — v and +v where 
is the optical depth at frequency v from the centre to the 
outer edge of the ejecta. The line profile is therefore given 
by 


Tv) 

I(-v) 


exp(—2r^v) 


( 11 ) 


(i) i{u) ~ u 7 — 1 (a > 0, 2/3 — 3 + a > 0) 


Lucy et al. (1989) presented several examples for both 
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Figure 5. Set of models with i(r ) oc r 4 (i.e. /3 = 2.0), R in /Rout — 0.2, v(r) oc r and t; max — 1 illustrating the effects of varying the 
dust optical depth r and albedo lj. Peak fluxes are scaled to unity. The black dotted line marks v = 0 and the red dotted lines mark 
^min and +n min . 


the analytical case of the perfect absorber and a Monte Carlo 
model for grains with u> = 0.6. We present the same cases in 
Fig. 4 and note that the resulting profiles exhibit the same 
features and shape. Of particular interest is the scattering 
wing that appears beyond the maximum velocity (i> m ax = 1) 
on the red side of profiles in the partial scatterer case as a 
result of the packets doing work on the expanding sphere. 
This was noted by Lucy et al. (1989) as a potential diagnos¬ 
tic for the presence of dust in the ejecta of a supernova and 
we will discuss this further in Section 4.3. 

4.3 The sensitivity of the variable parameters 

It is of general interest to establish potential diagnostic sig¬ 
natures in the line profiles of supernovae and their rem¬ 
nants in order to trace dust formation more effectively. We 
here discuss the effects of the main parameters of interest, 
namely: 

• the maximum velocity Vmax 

• the ejecta radius ratio -Rm/I?out 


• the dust optical depth r 

• the dust albedo ai 

• the density profile index /3, where p oc r 


4-3.1 The maximum velocity Vinax 

The maximum velocity is defined as the velocity at the outer 
edges of the line emitting region for a given line. The max¬ 
imum velocity may vary between different spectral lines or 
doublets due to different locations of species having differ¬ 
ing ionization thresholds. Clearly, the larger the maximum 
velocity used the wider the profile becomes. To some extent 
therefore the steepness of the density profile and the maxi¬ 
mum velocity can act to counter each other since a steeper 
density distribution narrows the profile (see Section 4.3.5). 
The shape of the wings of the profiles, however, generally 
precludes much degeneracy in this aspect - the overall shape 
of the line profile can be used to determine the exponent of 
the density distribution to within a relatively small range. 

More important is the effect that the maximum veloc- 
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Figure 6 . Set of models with i(r ) oc r ~for 0 = 1.0 (left), 0 = 1.5 (middle) or 0 = 2.0 (right), ui = 0, Ri u /R 0 ut = 0.2, v (r) oc r and 
n m ax = 1 illustrating the effects of varying the dust optical depth r. Peak fluxes are scaled to unity. The black dotted line marks v = 0 
and the red dotted lines mark —u m ; n and +?; Tn j T1 . 


ity has on the overall optical depth. Since the outer ra¬ 
dius is calculated directly from the maximum velocity (as 
Rout = Vmax x t where Vjnax is determined from the blue 
side of the observed line profile), the overall volume of the 
ejecta is determined solely by this value and the ratio of the 
inner and outer radii. The total dust optical depth to which 
the radiation is exposed can therefore be greatly affected by 
even a relatively small change in the maximum velocity for 
fixed values of the other parameters. Practically, however, 
the maximum velocity can usually be fairly well determined 
from the observations (identified as the point where the flux 
vanishes on the blue side) and may be further constrained 
through modelling. 


4-3.2 The ejecta radius ratio Ri n /R out 

As already discussed in Section 4.1, the width of the flat top 
is determined by the ratio of the inner and outer radii, the 
exponent of the velocity profile and the maximum velocity. 
We assume that the supernova is in free expansion from 


just a few months after the explosion and therefore r = vt 
such that within the ejecta the velocity profile takes the form 
v oc r at a fixed time i.e. the supernova expands self-similarly 
(Li & McCray 1992; Xu et al. 1992; Kozma & Fransson 
1998b). For this case, Riu/Rout is given by 


-rtout v max 

where it is often possible to constrain Knm and Vjnax to a 
relatively narrow range simply from the observed line profile. 

The majority of spectral lines emitted from supernovae 
and supernova remnants are expected to have a flat top be¬ 
fore dust attenuation effects since it is rare for these objects 
to form a completely filled nebula. However, even a very 
small amount of dust attenuation may result in the line pro¬ 
file appearing to be smoothed at its peak. 

The effects of absorption by dust on a line profile for a 
filled nebula with Ri n /Rout = 0, as opposed to a detached 
shell, are shown in Fig. 4. All profiles have been scaled to 
unit flux at their peaks. 
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Figure 7. The variation of albedo and extinction efficiency (Q e xt) with grain radius at A = 656 nm for Zubko et al. (1996) BE amorphous 
carbon, Draine &; Lee (1984) astronomical silicate and the MgSiOa and MgFeSiCH samples of Jager et al. (2003) and Dorschner et al. 
(1995), respectively. A linear grain size scale is presented on the left and a log scale on the right. 


^.3.3 The dust optical depth t (detached shell case) 

As expected, greater attenuation of the original line profile 
is seen on the red side (see Figs 5 and 6 ). The profiles are 
most revealing at lower dust optical depths since the effects 
of the asymmetric absorption can be seen in different sec¬ 
tions of the profiles and the profiles therefore tend to exhibit 
more features. The region of the profile that is most clearly 
affected by dust absorption is the flat-topped region. A small 
amount of absorption in this region results in a skewed pro¬ 
file, with a fraction of the flat-topped section removed. The 
peak becomes blueshifted as a result, but only to the original 
value of — Knm, the minimum velocity corresponding to i?i n . 
In addition to the attenuation in this region, the red wing 
of the profile is also somewhat reduced, and the blue wing 
somewhat increased relative to their original symmetric po¬ 
sitions. The result is a relatively “jagged” looking profile, 
often with sharp changes at ±Vjnin- The profile is generally 
asymmetric, although the degree of absorption in the flat- 
topped region may sometimes make it seem as though the 
profile is in fact symmetric and uniformly blueshifted (see 
Section 4.5 for further discussion). Observationally, these 
sharp features might become smoothed due to insufficient 
spectral resolution. 


At high dust optical depths or when the ratio of the 
inner and outer radii is small, the entire profile is shifted 
to the blue and the peak moves beyond — Vjnin further into 
the blue. The profiles also tend to become more smooth and 
featureless. A set of models showing the effects of varying 
optical depths for different density profiles and dust albedos 
are presented in Figs 5 and 6 with -Rm/.Rout =0.2. 


4-3-4 The dust albedo lj 

In the past, there has often been a focus on the effects of 
absorption by dust on the shapes of line profiles and less 
attention has been paid to the potential effects of scattering 
by dust. In fact, line profiles can be significantly affected by 
scattering of radiation. The greater attenuation of radiation 
received from the receding portion of the ejecta results in an 
asymmetry of the line profile whereby the majority of the 
observed emission is located bluewards of the peak. How¬ 
ever, the effects of repeated dust scattering events within 
the ejecta can substantially alter the shape of a line profile 
and potentially can act to counter the blueshifted asymme¬ 
try. 

Not only does repeated scattering of photons increase 
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the number of potential opportunities for a given photon 
to be absorbed but it also results in continuous shifting of 
the frequency of the photon to the red. The photon must 
do work on the expanding shell of dust in order to escape 
and thus many of the photons are reprocessed beyond the 
theoretical maximum velocity on the red side of the profile. 
Even in the case of dust grains with a relatively low albedo, 
a surprisingly persistent wing on the red side of the profile is 
seen, generally beyond the maximum theoretical velocity of 
the emitting region. In the case of strong dust scattering and 
high dust optical depths, this can actively result in a shift 
in the overall asymmetry of the profile, with the majority of 
the emission being emitted redwards of the peak. The peak 
however, remains blueshifted (for example the bottom left 
panel of Fig. 5) or central (for example the bottom right 
panel of Fig. 5). For the line prohle to exhibit this effect 
requires the dust to be a nearly perfect scatterer; of the 
albedos plotted in Fig. 7 only the nearly transparent MgSiOs 
sample of Jager et al. (2003) exhibits such a behaviour. 

The combination of relatively low dust optical depths, 
initially flat-topped profiles, greater attenuation on the blue 
side along with increased flux on the red side due to scatter¬ 
ing can result in a profile that sometimes ends up appearing 
almost symmetrical, particularly if contaminants, such as 
narrow lines or blending with other broad lines, are present 
or if the resolution of the data is low. The potential for 
apparently symmetrical profiles that appear to have been 
uniformly blueshifted should be noted (see Figs 5 and 6 for 
examples of this). 

4-3.5 Density profile p oc r^ 13 

Whilst the density prohle of the dust may have some effect 
on the resulting profiles, it is the initial emissivity prohle 
(dependent on the gas density prohle) that has the greatest 
effect on the shape of the line prohle. In general, the steeper 
the emissivity distribution, the narrower the line prohle be¬ 
comes. The sides of the line prohle may become almost ver¬ 
tical for a very steep distribution since the majority of the 
emission then comes from a very narrow velocity range (see 
Figure 3). 

The dependence of the shape of the line prohle on the 
emissivity distribution is described analytically in Section 
4.1 for the case of very optically thin dust. However, for 
even fairly low dust optical depths, the density prohle plays 
a significant role in determining the shape of the line pro¬ 
file where it is affected by dust absorption. As previously 
discussed, at relatively small optical depths for reasonable 
f?in/7?out, a section of the hat-topped region is removed re¬ 
sulting in a peak at —Rmin. The shape of the prohle in this 
region is signihcantly affected by the density prohle. Shallow 
density prohles (low fi) produce a virtually linear variation 
in flux between — V m in and +Knm (for example the prohles 
in the left-hand column of Fig. 6). For a hxed dust opti¬ 
cal depth, the steeper the distribution becomes, the more 
concave the prohle becomes between — Vmin and -l-Vmin, ul¬ 
timately resulting in a clear shoulder to the prohle at +Vmin 
(for example the prohles in the right-hand column of Fig. 6). 
For extremely steep density distributions this can result in a 
double peaked prohle with a trough to the red of V = 0. An 
illustration of the effects on the line prohles of varying ft and 
r is shown in Fig. 6. As previously noted, these features may 


not be apparent in observed line prohles with poor spectral 
resolution. 


4.4 Inferring properties of the dust from the 
models 

The presence of an extended red wing at large positive ve¬ 
locities in combination with increased extinction on the red 
side at smaller positive velocities can allow the values of r 
and u to be well constrained. Where this occurs it is possible 
to translate these values into a dust mass and grain size for 
a given species or combination of species using grain optical 
properties and Mie theory (see Figure 7). 

For amorphous carbon, the albedo generally increases 
with grain size. The presence and extent of any scattering 
wing on the red side of the observed prohle can therefore 
help to place limits on the grain radius. However, the greater 
the grain radius used the smaller the available cross-section 
for interaction per unit dust mass. Larger masses of dust 
are therefore required to ht the same degree of absorption 
if a larger grain size is used. This is in contrast to SED ra¬ 
diative transfer modelling where larger grain sizes generally 
result in less dust being required to ht the IR portion of the 
SED (W15). These two techniques in tandem may therefore 
provide limits on grain sizes for different species or combi¬ 
nations thereof. 

It is known that the use of different optical properties 
may substantially alter dust masses derived using SED fit¬ 
ting for a given species of specihc grain size (e.g. Owen & 
Barlow (2015)). However, the use of different sets of grain 
optical constants in our models seems to have only a minor 
effect on the required dust masses, except for cases where 
the albedo is close to unity (pure scattering grains). 

4.5 The wavelength dependence of dust 
absorption 

The greater the dust optical depth, the more attenuation 
of the line there will be. As expected, the red side of the 
prohle suffers a greater degree of absorption than the blue 
side. The resulting asymmetry is somewhat more complex 
than perhaps previously thought. Dust has repeatedly been 
cited as the agent responsible for the apparent blueshifting 
of supernova line prohles in the manner of the prohles pre¬ 
sented in Fig. 4; that is, relatively high optical depths result 
in an overall shift of the entire prohle towards the blue. The 
relationship between the blueshifting of the peaks of prohles 
and their wavelength has been discussed by several authors 
in relation to dust formation (Smith et al. 2012; Fransson 
et al. 2014; Gall et al. 2014). 

In practice, a relatively large dust optical depth is re¬ 
quired to actively shift the peak of the prohle bluewards of 
its natural — V m i n position (corresponding to the velocity at 
the inner radius of the shell) unless this value is very small 
in comparison to Knax be. the prohle originally had a very 
narrow hat top. In many cases, it seems likely that the dust 
may not be optically thick and the blueshifting of the peak 
of the prohle is just a result of attenuation in the Hat-topped 
section (close to i?i n ). The peak would then tend to be lo¬ 
cated at — Hmin- 

Since dust absorption is wavelength dependent for 
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Figure 8. Model line profiles for Ha (6563A in red), H/3 (486lA in yellow) and Pa<5 (10049A in blue) for optically thin (upper) and 
optically thick (lower) cases respectively. All models adopted density profile p(r) oc r~ 4 (i.e. /3 = 2), velocity profiles v(r) oc r and radii 
ratio -Ri n /-R ou t = 0.2. The grain radii used were a = 0.001 pm (left), a = 0.1 pm (middle) and a = 1.0 pm (right). All the above models 
used Zubko et al. (1996) BE amorphous carbon. 
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Figure 9. The variation of amorphous carbon dust absorption efficiency with grain size. The grain radii plotted are a = 0.001 pm (left), 
a = 0.1 pm (middle) and a = 1.0 pm (right). The vertical lines mark the wavelengths of Ha (6563A in red), H/3 (486lA in yellow) and 
Pa<5 (10049A in blue). 


2-ko, < A, one might expect the position of the peak line 
flux to be dependent on the wavelength of the line being 
considered. We note here that whilst variations of the peak 
velocity of a line as a function of line wavelength may oc¬ 
cur in cases of high dust optical depths or small Rm/Rout, 
this may not be the case for many supernova lines emitted 
from ejecta with low dust optical depths. The wavelength- 
dependence of dust absorption instead can result in differing 
degrees of extinction in the flat-topped region of each profile 
but still leave the peak at its blueshifted position of — Vmin- 
If this is the case then there would be no reason to expect 


a variation in the position of the peaks of profiles to be cor¬ 
related with the wavelength dependence of dust absorption. 
Instead one would expect it potentially to trace the location 
of different ions within the ejecta, possibly with different 
Vmin values observed for different species. 

For lines from the same ion, for example the Balmer and 
Paschen lines of Hi, we might expect to see peaks at the same 
position but differing degrees of absorption. At high spectral 
resolutions, it might be possible to detect differences in the 
shapes of the line profiles, particularly between —Vmin and 
TVmin where the steepness of the incline traces the degree 
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velocity (10 4 km s 4 ) 


Figure 10. Amorphous carbon smooth dust fit to the day 714 
Ha line of SN 1987A using an MRN size distribution, illustrating 
the underestimation of the red scattering wing for small grain 
radii. Model parameters are the same as the smooth dust fit for 
day 714 (Table 4) except for the grain size distribution and dust 
mass: M t j uat = 8.0 X 10 _6 Mq, a m ; n = 0.005 /mi, a max = 0.25 /tm 
and n(a) oc a -3 ' 5 . 

of dust absorption. This can be seen in Figure 8 where we 
illustrate the effects of the wavelength dependence of dust 
absorption for three lines, Ha (6563A), H/3 (4861A) and 
Pad (10049A). All lines were modelled using three different 
grain sizes and for both optically thin and thick dust cases. 
We also show the variation of the absorption efficiency with 
wavelength for three different amorphous carbon grain sizes 
in Fig. 9. 


5 RESULTS FOR SN 1987A 

We have modelled the Ha line of SN 1987A at days 714, 806, 
1862, 2211, 2875, 3500 and 3604, and the [O 1 ] A6300,6363 A 
doublet at days 714, 806, 1054 and 1478. After day 3604 the 
Ha profile begins to become dominated by emission from the 
reverse shock and the structure of the emitting region may 
no longer be approximated by a single shell model as we do 
here (Fransson et al. 2013). The [O 1 ] AA6300,6363 A doublet 
becomes too weak to model after day 1478 (see Fig. 1). We 
continue to adopt a velocity profile V(r) = j^ max r and treat 
the variable parameters listed at the start of Section 4.3. 
Whilst the albedo and optical depth are not varied directly, 
they are altered by adjusting the dust mass, Afdust, and 
the grain size, a, which together determine the albedo and 
optical depth via Mie theory and the optical properties of 
the dust. 

I 11 all models, the ejecta occupies a shell with inner ra¬ 
dius -Rm and outer radius R ou t- Packets are emitted accord¬ 
ing to a smooth density profile assuming recombination or 
collisional excitation such that i(r) oc p(r) 2 oc r~' 2,i . Initially 
the dust is considered to have a smooth density distribution 
and is assumed to be coupled to the gas so as to follow the 
same radial profile. A clumped distribution of dust is con¬ 
sidered later (see Section 5.2). 


Table 3. Observed luminosities of the Ha line and estimated 
electron scattering optical depths from iij n to 7? out for the radii 
detailed in Tables 4 and 5 based on an assumed gas temperature 
of 10,000 K. 


Ha 

day L oba 

(10 37 erg s _1 ) 

-^-'undep/ 

-^'obs 

[Oi] 

-^obs 

(10 37 erg s -1 ) 

-^undep/ 

-^obs 

T e 

( 10 - 2 ) 

714 

1.36 

1.65 

0.313 

3.57 

1.44 

806 

0.57 

1.77 

0.0942 

3.57 

0.840 

1054 



0.0242 

3.23 


1478 



0.00185 

2.70 


1862 

0.0063 

2.06 



0.159 

2211 

0.0041 

2.07 



0.0378 

2875 

0.0019 

2.84 



0.0219 

3500 

0.00079 

3.16 



0.0125 

3604 

0.00098 

3.27 



0.0149 


Assuming an electron temperature of 10,000 K, we es¬ 
timated the total electron scattering optical depths between 
Rin and R ou t based on the observed fluxes of the Ha recom¬ 
bination line. A temperature of 10,000 K for the recombining 
material is likely too high at the epochs considered but we 
adopt it in order not to underestimate electron scattering 
optical depths. The values we calculate from the observed 
Ha luminosities are listed in Table 3. Since the total elec¬ 
tron scattering optical depths at these epochs are negligibly 
small we therefore do not include electron scattering in the 
models. 

There is rarely a unique set of parameters that provide 
the best fit to the data. However, the majority of the param¬ 
eters of interest can be well constrained from our modelling 
by considering different elements of the shape of the profile. 
In particular, by constructing fits to the data using minimum 
and maximum limits for the grain radius, credible lower and 
upper bounds on the dust mass formed within the ejecta 
may be derived. We present here fits to the data obtained 
using both small and large values of the grain radius a since 
it is the grain size which has the most significant effect on 
the overall dust mass required to reproduce the line profile 
(see Section 4). 

All of our models are of a dusty medium composed 
solely of amorphous carbon grains. We use the optical con¬ 
stants from the BE sample presented by Zubko et al. (1996). 
Although previous SED modelling of SN 1987A limited the 
fraction of silicates present in the dusty ejecta to a maxi¬ 
mum of 15% (Ercolano et al. (2007), W15), the recent work 
of Dwek & Arendt (2015) has suggested that a large mass 
of mostly silicate dust may have formed at early epochs (~ 
615 d). It is therefore useful to consider the effects on our 
models of using silicate dust. We discuss this in detail in 
Sections 5.7 and 5.8. 

For each profile, the maximum velocity is initially iden¬ 
tified from the data as the point where the emission vanishes 
on the blue side and is then varied throughout the modelling 
in order to produce the best fit. The equivalent point on the 
red side is indeterminate from observations due to the effects 
of dust scattering. We determine the approximate value of 
Vmin by examining the width of the profile near its peak. 
Using the features and shapes presented in Figs 5 and 6 as 
a guide, we first examined the observed profile for any obvi- 
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velocity (10 4 km s 1 ) 





Figure 11. Best model fits to the SN 1987A Ha line at day 714 and day 806 for the parameters detailed in Tables 4 and 5. The two 
fits on the left are smooth dust models using amorphous carbon grains of radius a = 0.35 pm and the two fits on the right are clumped 
dust models using amorphous carbon grains of radius a = 0.6 pm. 


ous points of inflection or abrupt changes in the steepness of 
the profile. If these were observed then they were compared 
to similar changes in theoretical profiles which allowed us 
to estimate the value of Vmin■ If none were observed, then 
a model setting VAin to be the velocity of the profile peak 
was considered. Where neither of these approaches yielded 
a good model (this was rare) we iterated over a range of 
values of as with other variable parameters such as 

the dust mass. On the red side the theoretical minimum ve¬ 
locity often falls at a similar velocity to the 6583A line so 
any dust-induced features near this wavelength that would 
allow a more accurate determination of V m in can be over¬ 
whelmed by the nebular line. Having determined the min¬ 
imum and maximum velocities, the ratio of the inner and 
outer radii of the supernova ejecta can be determined since 
-Rin/-Rout = Vmin/hmax• The outer radius is calculated from 
the epoch and the maximum velocity. 

The only parameters that remain to be determined are 
the exponent of the density profile /3, the mean grain ra¬ 
dius and the total dust mass. The shape of the blue wing is 
solely a product of the density profile and the dust mass; the 
height and shape of the red wing is a product of these and 
also of the scattering efficiency of the grains (the albedo <n); 
the extent and shape of the asymmetry in the flat-topped 
portion of the profile is a function of only the total dust opti¬ 
cal depth determined by the dust mass and the grain radius. 
By iterating over these three parameters, an excellent fit to 
the data can usually be obtained. 

Models are produced in the same manner for the 
[O i] AA6300,6363 A doublet as for the single Ha line, with 
each component of the doublet being modelled indepen¬ 
dently and the resulting profiles added according to a speci¬ 
fied ratio. Although the theoretical intrinsic flux ratio is 3.1 
for optically thin emission (Storey & Zeippen 2000), the ac¬ 
tual ratio between the two components can be affected by 
self-absorption (Li & McCray 1992) and we therefore left it 
as a free parameter. The deduced doublet ratios are listed 
in Tables 4, 5 and 6. 

For all lines, though particularly at very late epochs, 
even small fluctuations in the adopted value of the contin¬ 
uum level can have a substantial effect on the fit to the 
resulting profile. Since it is not feasible to establish the level 
of the continuum so precisely, the value of the continuum 
has been left as a free parameter that may be adjusted (to 
within sensible margins) in order to allow for the widest pos¬ 
sible dust mass range to be determined. We generally find 


it is necessary to assume a continuum level that is slightly 
lower where the dust mass is higher. The [O l]A6300,6363 A 
doublets at days 1054 and 1478 are weak relative to the con¬ 
tinuum and are also blended with the wings of other lines 
making it difficult to fit their wings accurately. We aim to fit 
the lines between approximately —3000 km s _1 and +5000 
km s -1 but present a wider velocity range for context (for 
example see Fig. 13). 

Fits to the Ha line profile at days 2211 and 3500 are 
omitted for the sake of space but are very similar to those 
of days 1862 to 3604. All profiles have been smoothed to 
approximately the same resolution as the observed profiles 
using a moving-average procedure. Parameters for the mod¬ 
els at all epochs including days 2211 and 3500 are detailed 
in Tables 4 to 6. 

5.1 Smooth Density Models for SN 1987A 

Even at the earliest epochs there is a substantial wing on 
the red side of the Ha line profile that cannot be fitted 
by scattering from moving grains with a low albedo. The 
minimum required albedo is approximately u> ~ 0.5 imply¬ 
ing relatively large grain radii. As previously discussed, the 
larger the grain size the larger the mass of dust required to 
reproduce the same optical depth. Fig. 10 illustrates the fit 
for the day 714 Ha profile for the case where a classic MRN 
(Mathis et al. 1977) grain size distribution is adopted, with 
a m i n = 0.005/rm, a max = 0.25 fim and n(a) oc a _3 '°. It can 
be seen clearly that the extended red wing is significantly un¬ 
derestimated. Since the albedo of amorphous carbon grains 
varies significantly with grain radius (see Fig. 7) we can es¬ 
tablish a strong lower bound to the mean dust grain radius, 
which we estimate to be a > 0.35 pm. This is the smallest 
grain size that is still capable of reproducing the red scatter¬ 
ing wing at all epochs and we therefore use this lower limit 
value throughout our smooth density modelling. 

The inner and outer radii of the ejecta are calculated at 
each epoch from the maximum velocity used, the day num¬ 
ber and the specified ratio Ri n /Rout ■ The radii generated are 
consistent with those used in previous models of SN 1987A 
(Ercolano et al. (2007), W15) and the minimum velocities 
for both the [O i] and Ha line emitting regions are rela¬ 
tively consistent with those obtained by Kozma & Fransson 
(1998b) who estimate that hydrogen extends into the core 
to a depth of < 700 km s _1 and the oxygen reaches down 
to ~ 400 km s - They are also consistent with predictions 
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Figure 12. Best model fits to the SN 1987A Ha line at days 1862, 2875 and 3604 for the parameters detailed in Tables 4-6. On the 
top row are smooth model fits with amorphous carbon grains of radius a = 0.35 /im. On the middle and bottom rows are clumped model 
fits with amorphous carbon grains of radii a = 0.6 /im and a = 3.5 /im, respectively. 


from 3D explosion models at the time of shock-breakout 
that predict the oxygen to reach to a depth of ~ 200 km s _1 
(Hammer et al. 2010; Wongwathanarat et al. 2015). Figs 11 
to 13 show the best fits to the data for days 714 to 3604 
whilst Table 4 details the parameters used. 

It can be seen from Tables 4-6 that, in order to repro¬ 
duce the blueshifts seen in the [O i] AA6300,6363 A doublet, 
considerably larger dust masses are required than to fit the 
Ha line at the same epoch. Although the same maximum 
velocities and therefore outer radii are used in our [O i] and 
Ha models, the inner radii for the [O i] models are signif¬ 
icantly smaller and the density distribution much steeper. 
This implies that [O i] is concentrated towards the centre 
of the ejecta whereas Ha is more diffuse. This is broadly in 
agreement with 3D explosion dynamics models that suggest 
that a few hours after the explosion the heavier elements 
will, in comparison to hydrogen, be located more centrally 
in the ejecta with “bullets” of heavier material reaching the 
outer edges (Hammer et al. 2010). If dust is forming in the 
inner regions of the ejecta then the majority of the [O i] emis¬ 
sion must travel through the newly formed dust whereas the 
more diffuse Ha emission has a greater chance of escaping 
unaffected. This may explain the difference between the dust 
masses needed for the [O i] and Ha models. 


5.2 Clumped Dust Models for SN 1987A 

A number of investigators have presented arguments for the 
material in the ejecta of SN 1987A being clumped (Lucy 
et al. 1991; Li & McCray 1992; Kozma & Fransson 1998b) 
and so we consider clumped models for the ejecta dust to 
be more realistic than smoothly distributed dust models. It 
has been shown through the modelling of optical-IR SEDs 
that when dust is assumed to have a clumped distribution 
then the derived dust masses can be significantly larger than 
for the case of dust that is distributed smoothly between 
the inner and outer radii (e.g. Ercolano et al. (2007); Owen 
& Barlow (2015)). We present two sets of fits to the line 
profile based on the clumped dust modelling of W15, one 
set with a minimum grain size and one set with a maximum 
grain size. Each fit is based on the best fitting smooth model 
such that the photon packets are emitted assuming a smooth 
radial density profile. However, the dust is no longer coupled 
to the gas but instead is located entirely in clumps of size 
Rout/25. The clumps are distributed stochastically between 
Ri n and R ou t with the probability of a given grid cell being a 
clump proportional to r~^ where i(r) oc r~ 213 . The number 
of clumps used is determined by the clump filling factor / 
which is kept constant at / = 0.1. All properties are fixed 
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Table 4. The parameters used for the best fitting smooth models of SN 1987A with amorphous carbon grains of radius a = 0.35 fi m. 
Optical depths are given from R; n to R :m \. at A = 6563 A for Ha and A = 6300 A for [O i]. Values of Ty are very close to the quoted values 
of r Ha ■ 



day 

^max 

(km s -1 ) 

^min 
(km S -1 ) 

Rin/ Rout 

P 

Afdust 

(M @ ) 

-Rout 

( cm ) 

Rin 

(cm) 

[O i] ratio 

T\ 

[O I] 

714 

3250 

228 

0.07 

2.9 

9.65x 10 -5 

2.00X10 16 

1.40X10 15 

2.6 

3.60 

[O I] 

806 

4000 

240 

0.06 

2.4 

1.50x 10 -4 

2.79X10 16 

1.67X10 15 

2.3 

2.86 

[OI] 

1054 

4300 

215 

0.05 

2.1 

2.35x 10 -4 

3.92X10 16 

1.96X10 15 

2.7 

2.23 

[O I] 

1478 

4500 

180 

0.04 

1.7 

2.95x 10 -4 

5.75X10 16 

2.30X10 15 

3.0 

1.30 

Ha 

714 

3250 

813 

0.25 

1.2 

2.10x 10 -5 

2.00X10 16 

5-OlxlO 15 


0.61 

Ha 

806 

4000 

880 

0.22 

1.9 

3.80x 10 -5 

2.79X10 16 

6.13X10 15 


0.59 

Ha 

1862 

8500 

1275 

0.15 

1.9 

5.00x 10 -4 

1.37x 10 17 

2.05X10 16 


0.35 

Ha 

2211 

9000 

1260 

0.14 

1.9 

9.25x 10 -4 

1.72X10 17 

2.41X10 16 


0.42 

Ha 

2875 

9500 

1330 

0.14 

1.9 

1.50x 10 -3 

2.36X10 17 

3.30X10 16 


0.36 

Ha 

3500 

10000 

1400 

0.14 

1.9 

3.35x 10 -3 

3.02X10 17 

4.23X10 16 


0.49 

Ha 

3604 

10250 

1333 

0.13 

1.9 

4.20x 10 -3 

3.19X10 17 

4.15X10 16 


0.55 


Table 5. The parameters used for the best fitting clumped models of SN 1987A with amorphous carbon grains of radius a = 0.6 fin i. 
Optical depths are given from R\ tl to f? ou t at A = 6563 A for Ha and A = 6300 A for [O i]. Values of Ty are very close to the quoted values 
of r Ha ■ 



day 

^max 

(km s -1 ) 

^rnin 
(km S -1 ) 

Rin/ Rout 

P 

Afdust 

(M @ ) 

Rout 

( cm ) 

Rin 

(cm) 

[O i] ratio 

T\ 

[O i] 

714 

3250 

228 

0.07 

2.7 

2.00x 10 -4 

2.00X10 16 

1.40X10 15 

2.3 

3.84 

[OI] 

806 

4000 

240 

0.06 

2.3 

4.00x 10 -4 

2.79X10 16 

1.67X10 15 

2.0 

4.02 

[O I] 

1054 

4300 

215 

0.05 

2.3 

7.50x 10 -4 

3.92X10 16 

1.96X10 15 

2.3 

3.85 

[O I] 

1478 

4500 

180 

0.04 

2.0 

l.lOx 10 -3 

5.75X10 16 

2.30X10 15 

2.8 

2.65 

Ha 

714 

3250 

813 

0.25 

1.4 

5.50x 10 -5 

2.00X10 16 

5-OlxlO 15 


0.87 

Ha 

806 

4000 

880 

0.22 

1.8 

9.00x 10 -5 

2.79X10 16 

6.13X10 15 


0.76 

Ha 

1862 

8500 

1190 

0.14 

1.9 

1.20x 10 -3 

1.37x 10 17 

1.91X10 16 


0.46 

Ha 

2211 

9000 

1260 

0.14 

1.9 

3.00x 10 -3 

1.72X10 17 

2.41X10 16 


0.73 

Ha 

2875 

9500 

1140 

0.12 

2 

8.00x 10 -3 

2.36X10 17 

2.83X10 16 


1.05 

Ha 

3500 

10000 

1200 

0.12 

2 

1.35x 10 -2 

3.02X10 17 

3.63X10 16 


1.08 

Ha 

3604 

10250 

1230 

0.12 

2 

1.70x 10 -2 

3.19X10 17 

3.83X10 16 


1.22 


Table 6. The parameters used for the best fitting clumped models of SN 1987A with amorphous carbon grains of radius a = 3.5 fim. 
Optical depths are given from R; n to R (m \. at A = 6563 A for Ha and A = 6300 A for [O i]. Values of Ty are very close to the quoted values 

of T Ha . 



day 

^max 

(km s -1 ) 

^min 
(km S -1 ) 

Rin/ Rout 

P 

Mdust 

(Mq) 

Rout 

( cm ) 

Rin 

(cm) 

[O i] ratio 

T\ 

[O i] 

714 

3250 

228 

0.07 

2.9 

1.50x 10 -3 

2.00X10 16 

1.40X10 15 

2.3 

4.20 

[O i] 

806 

4000 

240 

0.06 

2.3 

2.70x 10 -3 

2.79X10 16 

1.67X10 15 

2.1 

3.95 

[O I] 

1054 

4300 

215 

0.05 

2.3 

5.50x 10 -3 

3.92X10 16 

1.96X10 15 

2.5 

4.12 

[O I] 

1478 

4500 

180 

0.04 

1.9 

8.00x 10 -3 

5.75X10 16 

2.30X10 15 

2.8 

2.81 

Ha 

1862 

8500 

1190 

0.14 

1.9 

l.OOx 10 -2 

1.37x 10 17 

1.91X10 16 


0.55 

Ha 

2211 

9000 

1260 

0.14 

1.9 

2.40x 10 -2 

1.72X10 17 

2.41X10 16 


0.85 

Ha 

2875 

9500 

1140 

0.12 

2 

6.00x 10 -2 

2.36X10 17 

2.83X10 16 


1.15 

Ha 

3500 

10000 

1200 

0.12 

2 

1.15x 10 _1 

3.02X10 17 

3.63X10 16 


1.34 

Ha 

3604 

10250 

1230 

0.12 

2 

1.25x 10 _1 

3.19X10 17 

3.83X10 16 


1.31 


from the smooth models with the exception of the grain 
radius, density profile exponent (/3) and the total dust mass. 

Models were again constructed using the smallest pos¬ 
sible grain radius (a=0.6 /im in the clumped case) in order 
to derive minimum dust masses for clumped distributions. 
By considering the extent of the red scattering wing, upper 
limits to the grain size were also derived with the purpose of 
limiting the maximum dust mass at each epoch. By steadily 


reducing the grain radius from an initial value of 5 /im (mo¬ 
tivated by the maximum possible grain size derived by W15 
for their day 8515 model), we produced a set of models with 
a maximum grain radius of a = 3.5 fim. 

The increase in grain size from the smooth case to the 
clumped case is necessary in order to have a slightly larger 
albedo. Grains of radius a = 0.35 /tm do not reproduce the 
red side of the profiles well for a clumped medium. This is 
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Figure 13. Best model fits to the SN 1987A [O i] A6300,6363 A doublet at days 714, 806, 1054 and 1478 for the parameters detailed in 
Tables 4-6. On the top row are smooth dust fits with amorphous carbon grains of radius a = 0.35/tm. On the middle and bottom rows 
are clumped dust fits with amorphous carbon grains of radii a = 0.6/tm and a = 3.5//m respectively. 


because when the dust is located in clumps the radiation 
is subject to less scattering as well as to less absorption. 
The reduction in scattering appears not to be compensated 
for by the increased dust mass and a larger grain radius is 
therefore required, particularly at day 714. 

For all but the Ha line at days 714 and 806 a similar fit 
could be obtained with either a grain radius of a = 0.6 /tm 
or a = 3.5 /im (see Figs 11 - 13). However, for Ha at days 
714 and 806 even a small change to the grain radius from 
0.6 /tm resulted in a significantly poorer fit, either over¬ 
estimating or under-estimating the red wing. We therefore 
conclude that the dust mass estimates produced for the Ha 
lines at days 714 and 806 for a grain radius of a = 0.6 /.tm 
are the best Ha-based estimates of the dust mass at this 
epoch. 

In our subsequent analyses, we adopt the values derived 
from our clumped models. Details of the parameters used 
are presented in Tables 5 and 6 and the fits are presented in 
Figures 11 and 12. 


5.3 Goodness of fit 

We detailed at the start of Section 5 the process by which pa¬ 
rameters were constrained in order to obtain good fits to the 
data. These fits were judged both by eye and by minimizing 
the MSE between the model and the observed data for each 
line profile. For those interested in the sensitivity of the fits 
to various parameters, in Tables 7 and 8 we detail the mean 


square error (MSE) for the Ha profile at days 714 and 2875 
for a range of dust masses and density profile exponents. All 
other parameters were kept fixed at their best-fitting values 
for the clumped models of Ha with a grain radius a = 0.6/tm 
as in Table 5. The MSE is calculated as 


1 

N 




,/mod,i ) 


(13) 


where N is the number of data points, /obs.i is the ob¬ 
served flux at the i th data point and / mo d,t is the modelled 
flux at the i th data point. The MSEs were calculated be¬ 
tween —5000 and +7000 km s -1 for the day 714 Ha profile 
and between —8000 and +8000 km s -1 for the day 2875 
Ha profile. Note that the MSEs should only be compared 
between models for a given observed line profile and not 
between different line profiles since each observation is asso¬ 
ciated with a different inherent error. 

For day 714, we find that increasing or decreasing the 
total dust mass by a factor of two with all other parameters 
fixed causes a substantial increase in the mean square error 
(by factors of 23 and 8.6 respectively) effectively ruling out 
these values. For day 2875 a similar variation is seen but 
with the MSE varying by factors of 1.4 and 3.0 for each 
case. The narrower range of MSEs at day 2875 compared to 
day 714 is due to a noisier profile which results in a greater 
allowed range of good fits. The sensitivity of the goodness 
of fit to the dust mass and density profile is similar for the 
other modelled epochs. 


MNRAS 000, 1-27 (2015) 

























































































18 Antonia Bevan and M. J. Barlow 


Table 7. MSEs illustrating the variation in goodness of fit for the Ho line profile for a range of dust masses with other parameters fixed 
at their best-fitting values for the clumped model with a = 0.6/rm as detailed in Table 5. The MSE is calculated between —5000 and 
+7000 km s' 1 for the day 714 Ho profile and between —8000 and +8000 km s“ 1 for the day 2875 Ho profile. A factor of zero represents 
the dust-free model. The best-fitting model is italicized. 


multiple of best-fitting mass 


0 

0.1 

0.5 

1.0 

2.0 

10 

Day 714 MSE (10 13 erg cm 2 s x ) 

0.167 

0.133 

0.043 

0.005 

0.115 

1.15 

Day 2875 MSE (10 -15 erg cm -2 s _1 ) 

0.0791 

0.0604 

0.0258 

0.0182 

0.0563 

0.288 


Table 8. Mean square errors illustrating the variation in goodness of fit for the Ho line profile for a range of density profiles with other 
parameters fixed at their best-fitting values for the clumped model with a = 0.6/rm as detailed in Table 5. The MSE is calculated between 
—5000 km s^ 1 and +7000 km s -1 for the day 714 Ho profile and between —8000 km s -1 and +8000 km s -1 for the day 2875 Ho profile. 
The best-fitting model is italicised. 


density profile exponent (ft) 



1.0 

1.2 

14 

1.6 

1.8 

Day 714 MSE (10 13 erg cm 2 s x ) 

0.0328 

0.0117 

0.005 

0.0184 

0.0410 



1.6 

1.8 

2.0 

2.2 

2.4 

Day 2875 MSE (10 15 erg cm 2 s 3 ) 

0.0282 

0.0205 

0.0182 

0.0193 

0.0255 


5.4 The effects of clumping 

As in the case of SED radiative transfer models, the 
dust masses required to reproduce the observations in the 
clumped scenario are considerably higher than for the 
smooth scenario. The dust masses differ between our smooth 
models for a = 0.35 /im and clumped models for a = 0.6/rm 
by a factor of approximately 3. The dust mass estimates are 
even larger when comparing clumped a = 0.6 pm models to 
clumped a = 3.5 pm models at later epochs. This does not 
take into account the increase in grain radius between the 
two cases however. This increase accounts for a reasonable 
fraction of this difference. We estimate the effects of clump¬ 
ing alone to increase the required dust mass by a factor of 
approximately 1.5-2.0 from the smooth case. 

5.5 More complex models 

Where blueshifted lines are observed in the spectra of 
CCSNe it is often the case that the Balmer lines of HI are 
less affected than the [O i] lines (Milisavljevic et al. 2012). 
This may be due to a difference in the location or distribu¬ 
tion of the emitting elements; if the neutral hydrogen was 
diffusely distributed throughout the envelope but the oxy¬ 
gen was co-located with the dust in the core and in clumps 
then this could result in [O i] emission undergoing greater 
attenuation than Ha. This geometry would be in line with 
previous models of SN 1987A that suggested that the dust¬ 
forming regions are likely to include those which are oxygen- 
rich (Kozma & Fransson 1998a). Clearly, any model of dust 
formation in the ejecta of a CCSN must consistently repro¬ 
duce all of the line profiles at a given epoch. The models 
presented in this paper thus far have coupled the gas and 
dust distributions for a fixed clump volume filling factor and 
clump size. The Ha and [O i] models therefore require dif¬ 


ferent dust masses with the [O i] models usually requiring a 
dust mass ~ 4 times larger than the Ha models. 

We now present a model that reconciles this difference 
by additionally varying the clump filling factor, clump size 
and emissivity distribution. We assume that neutral hydro¬ 
gen is likely diffuse throughout the ejecta and so maintains a 
smoothly distributed power-law emissivity distribution be¬ 
tween R[ n and R ou t for Ha. However, we now assume that 
dust mostly forms in dense regions of high metallicity and so 
restrict the [O l]A6300,6363 A emission to originate largely 
from the dusty clumps with only a small fraction emitted 
from the inter-clump medium. As previously discussed, the 
greater the covering factor of the dust the greater the albedo 
required in order to reproduce the Ha red scattering wing. 
In order to obtain both the strong blueshifting of the [O i] 
line and the extended red scattering wing observed in Ha 
a small number of dense clumps were required along with a 
small mass of diffusely distributed highly scattering dust in 
the inter-clump medium. 

In order to fit both line profiles simultaneously, we re¬ 
quired a very high albedo (w > 0.8) that demanded the 
inclusion of some fraction of silicate dust. Amorphous car¬ 
bon grains alone are incapable of producing this level of 
scattering for any grain size. We adopted a grain radius of 
a = 0.6^im, the same as that used in our initial clumped 
models and we varied the relative proportions of amorphous 
carbon and MgSiC >3 in order to achieve the necessary albedo. 
The adopted grain densities were p c = 1.85 g cm~ 3 for amor¬ 
phous carbon grains and p s = 2.71 g cm~ 3 for MgSiOa. 
The resulting dust model for day 714 used 75% MgSiC >3 and 
25% amorphous carbon by cross-sectional area with a vol¬ 
ume filling factor fv = 0.1 and a clump size R a ut/5- 90% 
of the dust mass was located in clumps with the remaining 
10 % distributed smoothly between Ri n and R ou t according 
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Figure 14. Fits to the Ha and [O l]A6300,6363 A lines at day 
714 using the more complex dust model described in Section 5.5 
with a dust mass of 2.3 X 10 — 4 Mq. 


to a power law p oc r. Clumps were distributed stochas¬ 
tically with probability oc r -8 compared to r~ 2 '' in our 
standard models discussed earlier. Equal numbers of [O i] 
packets were emitted from each clump. The increased steep¬ 
ness of the density profile is required to compensate for the 
clumped packet emission relative to the previous smooth 
distribution. Since the clumps are distributed stochastically 
according to the density profile, less flux is emitted from 
the central regions in a clumped emission model than in a 
smooth distribution model (since there are gaps between the 
clumps). In order to obtain a sufficiently steeply rising line 
profile, the density profile must therefore be steepened in 
clumped emission models. The adopted value of f3 does not 
significantly affect the best-fitting values of the other pa¬ 
rameters of interest however. Ha was distributed smoothly 
according to a density power law p(r) oc r~ 3 . R ou t was the 
same for all components (i.e. clumped dust, diffuse dust, 
[O i] emission and Ha emission) and was calculated using 
a maximum velocity of 3250 km s“ . The inner radius was 
Rin = 0.07-Rout for all components except the smooth Ha 
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emission which was emitted between R ln = 0.25R ou t and 

Rout - 

The total dust mass used was Md us t = 2.3 x 10 -4 Mg. 
This dust mass is very similar to that derived from our orig¬ 
inal clumped models of [O i] using amorphous carbon grains 
of radius a = 0.6/xm. The slight increase over our amor¬ 
phous carbon dust mass of 1.5 x 10~ 4 M@ is largely due 
to the higher grain density of MgSiOs. At this grain radius 
amorphous carbon and MgSiC >3 have similar extinction effi¬ 
ciencies and so the change in species and geometry does not 
substantially alter the dust mass. We therefore adopt the 
[O i] dust masses in our further analyses and consider the 
differences in our derived dust masses between Ha and [O i] 
to be the result of the clumped emission of [O i]. 

Fits to both the [O l]A6300,6363 A and Ha lines for day 
714 using these parameters are presented in Fig. 14. 


5.6 The effect of a grain size distribution 

It is important to consider the potential effect on the dust 
mass of modelling a grain size distribution instead of a single 
grain size. For a grain size distribution the overall extinction 
cross-section, Cext, at a given wavelength is 

/ a max 

Qext{a)n(a)na 2 da (14) 

min 

where Q ex t (a) is the extinction efficiency for a grain size 
a and n(a) is the number of grains with size a. The overall 
extinction efficiency is then 


Qext - ( 15 ) 

The scattering cross-section Q aca is similarly calculated. 
As a result of these calculations, there is rarely a single grain 
size that has the same albedo and extinction efficiency as a 
size distribution. Modelling a size distribution may there¬ 
fore alter the deduced dust mass. Since the models are only 
sensitive to the overall optical depth and albedo, it is not 
possible to deduce the grain size range or distribution and 
only single grain sizes are investigated (as presented above). 

Whilst this apparently limits the scope of the results, 
it is useful to consider the extent to which different grain 
size distributions would alter the derived dust masses. By 
considering a number of grain radius ranges and adopting 
a power law distribution with a variable exponent, we may 
gain some insight into the effects of adopting a distribution 
rather than a single size. As discussed in Section 5.1, for a 
classical MRN power law (n(a) oc a~ 3 ' 5 ) with a wide grain 
radius range (a min = 0.001 pm to a max = 4.0 pm) the de¬ 
rived albedo is much too small to reproduce the required 
wing seen at early epochs. We therefore adopt an approach 
whereby, for a number of grain size ranges, we adjust the 
exponent of the distribution until the overall albedo is the 
same as that seen for the best fitting single grain radius 
for the clumped distributions. We may then approximately 
calculate the required dust mass as 


M d = Ms Q^Aas) x 


P max n(a)a 3 da 

Jfl min V ’ _ 

Qext (a)n(a)a 2 da 


(16) 
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Table 9. Dust masses for day 714 clumped models of the Ha 
line using different grain size distributions and 100% amorphous 
carbon. The final column shows the factor of increase over the 
dust mass for the single size model (M = 7 X 10 -5 Mg with 
a = 0.6 /im) and p is the exponent of the grain size distribution 
n(a) oc a~ p . 


Table 10. Dust mass conversion factors for single size models us¬ 
ing grains of 100% Zubko BE amorphous carbon or 100% Draine 
& Lee silicate at A ~ 656 nm. / is the factor by which the dust 
mass changes on going from amorphous carbon to silicates. 


^min 

(/im) 

&max 

(/tm) 

P 

M 

(M q ) 

M/AIo.6 

0.001 

4.0 

2.45 

1.93 xlO -4 

2.76 

0.01 

4.0 

2.45 

1.93 xlO -4 

2.76 

0.05 

4.0 

2.52 

1.84 xlO -4 

2.62 

0.1 

4.0 

2.72 

1.61 xlO -4 

2.3 

0.5 

4.0 

8.20 

7.23 xlO -5 

1.03 


where the subscript s represents the single grain size 
quantities and the d subscript represents quantities for the 
grain size distribution. 

We calculate the required dust masses for the clumped 


a 

(/tm) 

carbon 

CO 

Qex.t 

a 

(/tm) 

silicates 

CO 

Qext 

Msil/Mamc 

0.6 

0.56 

2.61 

0.0583 

0.58 

0.08 

5.37 

0.6 

0.56 

2.61 

4.00 

0.56 

2.18 

13.0 

3.5 

0.62 

2.21 

0.0641 

0.64 

0.10 

0.65 

3.5 

0.62 

2.21 

1.020 

0.63 

2.15 

0.49 

3.5 

0.62 

2.21 

1.376 

0.62 

2.35 

0.61 


can calculate the mass of DL silicate that gives a fit equiv¬ 
alent to that for a single carbon grain radius. We consider 
the albedo for the grain radius needed for the best-fitting 


Ha model on day 714 for a selection of distributions with 
varying a m j n . These are presented in Table 9. It can be seen 
that in all cases, a larger dust mass is required for grain 
size distributions in order to reproduce the same profile as 
a single grain size. The conversion factors presented in the 
table are valid for any model with grain size a = 0.6 /im and 
may therefore also be applied to the models for day 806. 
We repeated the process for a = 3.5 /im but found that, 
in order to reproduce the required albedo, the distribution 
had to be heavily weighted towards the larger grains and 
that the value of a m i n had no effect on the required dust 
mass. Increasing the value of a m i n to larger values (> 2 /jm) 
does not have a significant effect either. This is because both 
extinction efficiency and albedo tend to a constant value 
with increasing grain radius and the adoption of different 
grain size ranges and distributions above a certain threshold 
results in only insignificant variations in these quantities. 

We conclude that if a distribution of grain sizes is indeed 
present, the deduced single size dust masses are likely to 
under-estimate the true mass of newly formed dust. 


amorphous carbon model, calculate the equivalent grain ra¬ 
dius for DL silicate that gives the same albedo and then 
calculate a new dust mass by allowing for the change in the 
extinction cross-section: 


/ Qamc \ 

/ &sil \ 

/ Psil \ 

' Qsil ' 

' flame ' 

V PamC > 


Because of the nature of the variation of albedo with 
grain radius for the Draine & Lee (1984) astronomical sili¬ 
cate (see Fig. 7), there is often more than one silicate grain 
radius that will give rise to the same albedo at a given wave¬ 
length. Some of the possibilities and the resulting mass con¬ 
version factors are given in Table 10. For our best fitting 
amorphous carbon models with a = 0.6 /tm (the first two 
entries in Table 10), using any fraction of silicates with ei¬ 
ther a = 0.6 /im or a = 3.5 /im would increase the dust mass. 
However, for the case of an amorphous carbon grain radius 
of a = 3.5 /an (the last three entries), using silicate dust 
would reduce the dust mass by a factor of about 2 relative 
to our amorphous carbon values. 


5.7 The effect of different grain species 

In our analyses so far, we have mostly focused on amorphous 
carbon as the species of interest. This was motivated by pre¬ 
viously published early epoch optical and IR SED analyses 
that found that the silicate mass fraction must be limited to 
<15% (Ercolano et al. (2007), W15). The recent suggestion 
by Dwek & Arendt (2015) that large masses of the glassy sil¬ 
icate MgSiCU may have formed at early epochs is discussed 
further in the next subsection. The parameters that affect 
the quantity of dust required by our models are the mean 
albedo and optical depth of the dust. There could be multi¬ 
ple combinations of grain species and sizes that result in a 
good fit to the data. 

We can evaluate the required change in dust mass when 
a medium of 100% silicates is used instead of amorphous 
carbon. Using the astronomical silicate optical constants of 
Draine & Lee (1984), which are ‘dirtier’ (with lower albedos) 
than the glassy pure MgSiOs sample of Jager et al. (2003). In 
a similar manner to the approach detailed in Section 5.6, we 


5.8 Modelling large masses of dust at early 

epochs: comparison with the results of Dwek 
& Arendt (2015) 

In a recent analysis of infrared SED data, DA15 suggested 
that it may be possible for a large mass (0.4Mg) of MgSiCU 
silicate dust to have been present in SN 1987A even at rel¬ 
atively early epochs (t ~ 615 d), since that species has very 
low IR emissivities. Up to this point, we have constructed 
models using Zubko et al. (1996) BE amorphous carbon dust 
but in the previous section we discussed the effect on derived 
dust masses of instead using Draine & Lee (1984) astronomi¬ 
cal silicate which has higher optical and IR emissivities than 
the glassy MgSiCU species considered by DA15. Our clump¬ 
ing structure in our models was based on that used by W15. 

We now consider models for day 714 based on the grain 
types used by DA15. We adopt a clumped structure equiv¬ 
alent to the preferred model of DA15 who considered 1000 
clumps with a filling factor of 0.09 and a negligible dust 
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Figure 15. Ha models using different grain species and dust masses. Models for the dust masses presented by Dwek &; Arendt (2015) 
are on the top and models using our minimum required dust masses are on the bottom. From left to right the dust species are composite 
grains (82% MgSiOa and 18% amorphous carbon by volume), pure MgSiOa, pure amorphous carbon and pure MgFeSiCU. A density 
distribution with f) = 2.3 was adopted with a filling factor / = 0.09 and an effective clump radius R e g/Rout = 0.044. All other parameters 
are the same as in Table 5. 


mass in the inter-clump medium. We calculate the effective 
spherical radius of our clumps by equating the volume of 
our cubic clumps to a sphere of radius R e g. Clumps of width 
■Rout/14 generate the desired R e a/Rout = 0.044 equivalent 
to that of DA15. In our code, using a filling factor of 0.09 
then generates 1034 clumps, similar to the number used by 
DA15. We ran a series of models (presented in Figs 15 and 
16) for both the Ha and [O l]A6300,6363 A line profiles. In 
each case we modelled the lines using a dust grain mixture 
as described by DA15 such that the medium comprised 18% 
amorphous carbon and 82% MgSiC >3 by volume. We adopted 
the same optical constants as used in their work (i.e. Jager 
et al. 2003 for MgSiCH grains and Zubko et al. 1996 for amor¬ 
phous carbon) and the same grain mass densities as DA15, 
p a = 3.2 g cm' 3 and p c = 1.8 g cm” 3 . In addition to mod¬ 
elling their composite grain case, we also considered three 
single species models, using Zubko BE amorphous carbon, 
MgSiC> 3 , and MgFeSiO .1 (in the latter two cases the optical 
constants were taken from Jager et al. 1994 and Dorschner 
et al. 1995). For each species we adopted the smallest single 
grain size that has an albedo of to « 0.6. The ejecta param¬ 
eters were as listed in Table 5, with the exception of the 
density distribution which we took to be p(r) oc r” 1 3 for 
Ha and p(r) oc r” 2 ' 3 for [O i] in order to optimise the best 
fits. 

For each species, two models are presented. The first 
adopts the minimum possible dust mass that provides a 
reasonable fit to the observed line profiles and the sec¬ 
ond uses the dust mass derived by DA15 for that specific 
species (M = 0.4 Mg for MgSiCH and M = 0.047 Mg for 
amorphous carbon giving a total composite dust mass of 
M = 0.447 Mq). We treated MgFeSiCH as we do the com¬ 
posite grains and adopted a dust mass of M = 0.447 Mq for 
it. Results from the models are presented in Figs 15 and 16. 

The [O i] models can display similar profiles for substan¬ 
tially different dust masses. This is a result of the relatively 


high optical depths within the clumps themselves. If a clump 
is optically thick then the majority of radiation that hits it 
will be absorbed and the profile becomes insensitive to how 
much dust is actually contained within the clump. For our 
[O i] minimum dust mass models, the optical depths within 
a clump over an effective clump radius R e g at 6300A are 
around r c i ump « 0.4. Over the entire nebula optical depths 
are very high and ~72% of the total flux is absorbed. Increas¬ 
ing the total dust mass therefore has only a small effect on 
the emergent line profile and once r c i ump > 1 then the line 
profile remains unchanged for increasingly large dust masses. 
It is because of this fact that we present only the smallest 
dust mass capable of reproducing the [O i] profiles seen in 
Fig. 16. The insensitivity of the [O i] profiles to dust mass is 
not the case for the Ha profile models (where r c i ump < 0.05 
for all of our models) and the Ha-fit dust masses presented 
in Fig. 15 therefore represent the most sensitive diagnostic 
of the dust mass for each grain type. All of our models dis¬ 
cussed in previous sections have significantly smaller clump 
optical depths (r c iu mp < 0.1), making them sensitive to dust 
mass variations. 

For all the [O i] line profile models, except for those 
using pure MgSiCH or pure Mg 2 Si 04 dust, the required dust 
masses are significantly less than those proposed by DA15. 
The [O i] profile obtained using DA15’s very large MgSiC >3 
dust mass of 0.4 Mq provides a reasonable fit, but the same 
dust mass significantly overestimates the blueshifting of the 
Ha line (Figure 15). We can place an upper limit on the 
mass of pure MgSiCH on day 714 of 0.07 Mq, as this is the 
highest mass for which a fit to the observed Ha profile can 
be obtained (Figure 16). 

Pure MgSiCH is extremely glassy, with very high albe¬ 
dos in the optical for a wide range of grain radii. At grain 
radii small enough to reduce the albedo to oj ~ 0.6, in 
order to fit the observed line profiles, the extinction effi¬ 
ciency in the optical becomes extremely low (see Fig. 7), 
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Figure 16. [O l]A6300,6363 A models using different grain species and dust masses. Models using the dust masses presented by DA15 
are on the top and models using our minimum required dust masses are on the bottom. From left to right the species are composite 
grains (82% MgSiOa and 18% amorphous carbon by volume), pure MgSiOa, pure amorphous carbon and pure MgFeSiO.i. A density 
distribution with /3 = 1.3 was adopted with a filling factor / = 0.09 and an effective clump radius R e ff/R ou t = 0.044. The ratio between 
the doublet components was 2.2. All other parameters are the same as in Table 5. 


with large masses of dust therefore required in order to pro¬ 
duce even a small amount of line absorption. However, for 
a given albedo, the extinction efficiencies increase by large 
factors if either carbon or iron is included in the grain. 
In the composite grain model the amorphous carbon com¬ 
ponent dominates the overall extinction due to its much 
larger extinction efficiency at small grain radii. Similarly, 
for MgFeSiCH (or Mgo.sFeo.sSiOs) grains the iron compo¬ 
nent leads to much larger optical and IR extinction efficien¬ 
cies and much lower dust mass upper limits. If the dust 
that formed at early epochs contained some fraction of ele¬ 
ments such as carbon, iron or aluminium, yielding ‘dirtier’ 
silicate grains or composite grains, then fits to the observed 
blueshifted line profiles imply low dust masses. We conclude 
that for dust masses as large as 0.07 M@ to have been present 
in SN 1987A’s ejecta as early as days 600-1000 then the dust 
would have to have been formed of glassy pure magnesium 
silicates. 


In order to be certain that there was no set of param¬ 
eters for which a dust mass of M = 0.447Mq comprising 
82% MgSiOs and 18% amorphous carbon by volume could 
result in a good fit, a thorough investigation of the variable 
parameters was performed. Having fixed the clump size, fill¬ 
ing factor, dust mass and composition as per the values de¬ 
tailed above and in DA15, we varied the density profile (/3) 
and grain radius a. Varying the maximum velocity and the 
ratio of the inner and outer radii was found to have little 
effect on the goodness of fit. The MSE for the Ha profile 
presented in the upper left panel of Fig. 15 was 0.599 (in 
units of 10 -13 erg cnh 2 s -1 ). This was improved to 0.246 
by increasing the grain radius to a = 0.6/rm and the density 
profile exponent to /3 = 1.5, which represents the best fit 
that we could achieve using the values described by DA15 
and a dust mass of M = 0.447Mq. However, the overall best 
fit we obtain for this scenario (see the lower left panel of 15) 


used a dust mass of M = 5 x 10 4 Mq giving a MSE=0.0058, 
substantially improving the fit. 


5.9 Unattenuated line fluxes 

The evolution of the SN 1987A Ha and [O l]A6300,6363 A 
line fluxes over time has been discussed previously by, for 
example, Li & McCray (1992), Xu et al. (1992) and Kozma 
& Fransson (1998b). We may use our clumped models to 
predict the unattenuated emitted line fluxes and consider 
their evolution through time. For each model, the fraction 
of the total line energy absorbed by the dust was predicted. 
We determined the total flux for each observed line profile 
and used the absorbed fraction from our clumped models 
for a = 3.5/rm to predict the undepleted flux of the line 
before attenuation by the dust. Gaps in the observed data 
due to contamination by narrow line emission were inter¬ 
polated over in order to estimate the flux of the broad line 
component. The observed Ha luminosities and predicted un¬ 
depleted luminosities are given in Table 3 along with the 
energy fraction absorbed by the dust in each model. No cor¬ 
rection has been made for interstellar extinction along the 
sightline to SN 1987A. There is very little change in these 
values if we adopt the models with a = 0.6 /im instead of 
a = 3.5 /j,m. Plots of the observed and undepleted line lumi¬ 
nosities are given for all modelled epochs of Ha and [O i] in 
Figure 17. 

We also present power-law fits to the time evolution 
of the unattenuated Ha and [O i] line fluxes. For Ha, we 
find that LHa(t) oc t -4 15 between days 714 and 3604. We 
can compare this value to the theoretical time dependence 
of the flux of a recombination line based on the dynam¬ 
ics of the ejecta for an environment in a Hubble-type flow 
r = vt. For a frozen-in ionization structure, the mean inten¬ 
sity of a recombination or collisionally-excited line per unit 
volume is locally proportional to the product of the densi- 
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ties of the recombining species i.e. Jhc, oc n e n p oc n 3 . The 
total luminosity of the line is therefore dependent on the 
volume V as Lhc, oc 1/V. Assuming a constant maximum 
expansion velocity, the luminosity should vary with time as 
L Ha (t) oc t~ 3 . 

This relationship is only true for a constant ionization 
fraction. This ‘freeze-out’ phase is estimated to have begun 
at ~ 800 d and first sets in at lower density high velocity re¬ 
gions, gradually moving inwards with time (Danziger et al. 
1991; Fransson & Kozma 1993). Since our modelling be¬ 
gins at day 714, the ionization fraction in the inner higher 
density regions is likely still decreasing due to recombina¬ 
tion during our first two epochs. This presumably accounts 
for the slightly steeper Lna(t) oc f -4 ' 15 that we find across 
all epochs. Kozma & Fransson (1998b) estimate that Ha 
emission from the outer regions begins to dominate over Ha 
emission from core regions for t > 900 d. If earlier epochs 
are ignored, the last five epochs (t > 1862 d) plotted in (Fig. 
17) exhibit a shallower trend that is in good agreement with 
the expected Lu a (t) oc t~ 3 evolution. 

The [O l]A6300,6363 A doublet exhibits a much steeper 
evolution, T[or] 0) oc t~ 1 ' 2 , than the Ha line (Fig. 17). These 
collisionally excited lines are very sensitive to the gas tem¬ 
perature, with emissivities that fall to low values for temper¬ 
atures below ~3000 K. The models of Li & McCray (1992) 
and Kozma & Fransson (1998a) predict that the gas tem¬ 
perature in the relevant [O i] emitting regions should have 
fallen below 1000 K after day ~1000. 


6 DISCUSSION 

Using Monte Carlo models that consider both the absorbing 
and scattering effects of dust, we have modelled the evolu¬ 
tion of the Ha and [O i] A6300,6363 A line profiles over time, 
enabling us to place constraints on the evolution of newly 
formed dust in the ejecta of SN 1987A. 

As can be seen in Fig. 12, even a small degree of asym¬ 
metry in observed supernova line profiles can be indicative 
of dust formation within the ejecta. In addition to this, a line 
profile that is consistently asymmetric through time requires 
increasingly large dust masses to account for a similar de¬ 
gree of blueshifting since the expansion of the ejecta would 
otherwise cause the dust optical depth to the edge of the 
ejecta to be reduced. 

In Section 5.8 we compared our results with those of 
Dwek & Arendt (2015) and concluded that large dust masses 
can only have been present at early epochs if the grains 
were formed purely of glassy magnesium silicates that con¬ 
tained no iron or carbon component and that even for pure 
magnesium silicates no more than 0.07 M 0 could have been 
present. We now compare our results with those of Lucy 
et al. (1989) and W15. 

Lucy et al. (1989) analysed the [O i] A6300,6363 A dou¬ 
blet for SN 1987A and estimated dust optical depths for a 
number of epochs. They translated these into dust masses 
for day 775 only. From our smooth flow modelling of the [O i] 
doublets, we obtain tv ~ 3.60 at day 714 and ry « 2.86 at 
day 806. These values are higher than the values given by 
Lucy et al. (1989) who derived tv = 1.19 at day 725 and 
tv = 1-25 at day 775. The value of the assumed albedo ac¬ 
counts for the majority of this discrepancy. Lucy et al. (1989) 



time (days) 



time (days) 

Figure 17. Predicted undepleted luminosities for the Ha line 
(above) and [O l]A6300,6363 A doublet (below) presented with 
the best power-law fit to the data. 

considered line profiles before and after dust condensation 
and concluded that any evidence of an extended red scat¬ 
tering wing was unconvincing. Accordingly, they adopted a 
model with perfectly absorbing dust (w = 0). For our amor¬ 
phous carbon models for the [O i] AA6300,6363 A profile 
using a grain radius a = 0.35/tm, we obtain an albedo of 
approximately uo = 0.5 at A = 6300 A. 

The dust masses derived by Lucy et al. (1989) at day 
775 (e.g. Mdust = 4.4 x 1O~ (, M 0 for amorphous carbon) are 
different to those obtained from our smooth dust modelling 
of the [O i] AA6300,6363 A doublet at day 806 (Md„ s t = 
1.5 x 1O _4 M 0 for amorphous carbon). There are three main 
reasons for the discrepancy. First, the albedo is significantly 
larger in our modelling as already discussed. A larger dust 
mass is therefore required to produce the same amount of 
absorption. Secondly, to match the extended red wing our 
required grain radius is considerably larger than the small 
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time (days) 

Figure 18. Derived dust masses for SN 1987A as a function of epoch. Red squares - dust masses derived by W15 from their photometric 
SED modelling of SN 1987A. Solid yellow line - W15’s sigmoid fit to their values. Dark and light blue asterisks - maximum (a = 3.5 //rn) 
and minimum (a = 0.6 /jm) dust masses respectively for the [O i] models for t < 1478 days and for the Ha models for t > 1862 days. 
Purple stars - predicted dust masses calculated as the mean of the maximum and minimum dust masses. Dashed green line - sigmoid fit 
to our predicted dust masses. 


grains (a < 0.1/im) adopted by Lucy et al. (1989). Larger 
grain radii reduce the total cross-section of interaction and 
so a greater dust mass must be present to compensate for 
this. Finally, the adopted maximum velocity (4000 km s -1 ) 
in our model is larger than the value adopted by Lucy et al. 
(1989, 1870 km s _1 ). The larger value of Vjnax increases the 
total volume of the ejecta significantly and therefore signif¬ 
icantly more dust is required to produce the same optical 
depth. 

Lucy et al. (1989) also noted that the dust optical depth 
increased rapidly after day 580 and that the rate of increase 
of the dust optical depth appeared to slow between day 670 
and day 775, the latest day that they considered. Our re¬ 
sults, for both clumped and smooth models, suggest that the 
dust optical depth actually drops between day 714 and day 
806 before starting to increase again at later epochs. This 
is consistent with the results of Lucy et al. (1989) where 
the slowing rate of increase of dust optical depth could be 
consistent with a turning point subsequent to day 775. 

We can also compare our dust masses with the mass 
estimates derived from SED-fitting by W15 (see Fig. 18). 


W15 used a sigmoid fit to their dust mass evolution, of the 
form 

M d (t) = ae be , (18) 

where a = 1.0 Mg (representing the limiting dust mass), 
b = —8.53 and c = —0.0004. Both their dust masses and 
this sigmoid fit are shown in Fig. 18. It exhibits an initial 
period of slow growth in mass followed by an intermediate 
period of accelerating growth followed by another slowing 
until a plateau is ultimately reached. In this sense it may 
be representative of the process of dust formation whereby 
initial conditions appropriate for grain growth gradually de¬ 
velop until optimal conditions are reached at an intermediate 
epoch when grain growth is at its fastest before conditions 
once again deteriorate and the rate slows again (as discussed 
by W15). Performing a least-squares regression to this func¬ 
tion using just our own derived clumped dust masses, we 
obtain a sigmoid fit with coefficients a = 1.0Mg, b = —10.0 
and c = —0.0004. These values are remarkably similar to 
those derived by W15. This sigmoid fit is also plotted in 
Fig. 18. 
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We find that at all epochs the dust masses derived by 
W15 are entirely within the dust mass ranges determined 
by our models. 

Our sigmoid fit to the mean of the maximum and mini¬ 
mum dust masses does not take into account any systematic 
effects of grain growth. At earlier epochs, whilst grains are 
still small relative to later epochs, the lower bound to the 
dust mass estimates may be more representative than the 
upper end; the reverse would be true at later epochs. This 
is in contrast to the sigmoid fit of W15, whose fits to their 
early epoch SEDs used an MRN distribution with grain radii 
between 0.005 and 0.25 fi m, whilst their fits to their last two 
epochs required grain radii between 3.005 and 3.25 /jm. The 
dust masses used for their sigmoid fit thus accounted for the 
effects of grain growth between the earlier and later epochs. 
As mentioned, we could not fit the extended red wings of 
the profiles at early epochs using an MRN distribution. W15 
found that at their earlier epochs they could not obtain SED 
fits with grain radii as large as ~ 1.0 /.an. However, they did 
not consider radii in between these size ranges, such as the 
grains with a « 0.6 /tm that we require at earlier epochs. For 
SED modelling it is generally the case that the larger the 
grain size used, the less dust is required to produce the same 
level of flux. This may account for the differences between 
W15’s earlier epoch dust masses and our own minimum dust 
mass estimates at similar epochs. The models of W15 used 
15% silicate dust, in contrast to our models which used 100% 
amorphous carbon dust. This could also contribute to the 
differences at early epochs, as could the use of different sets 
of optical constants - we used the BE amorphous carbon 
optical constants of Zubko et al. (1996) whereas W15 used 
AC constants from Hanner (1988). W15 found that in order 
to fit early epoch SEDs epochs (e.g. day 615) with Zubko 
ACH2 constants, smaller inner and outer ejecta radii were 
needed, with half as much dust (5.0 x 10~ 4 Mq) compared 
to the Hanner AC results. 

W15 derived a maximum possible grain size at late 
epochs, concluding that the grains could not be larger than 
~ 5 /im by day 8515. This is consistent with the maximum 
grain radii that we derive at our latest epochs. We find that 
grain radii most likely cannot have exceeded ~ 3.5 /im at 
day 3604 - the dust mass that we obtain using this grain 
radius is similar to the value predicted by W15’s sigmoid fit 
at that epoch. 

The relationship between ejecta dust grain radii and 
post-explosion time is important for understanding the like¬ 
lihood of dust surviving the passage of a reverse shock prop¬ 
agating back through the ejecta. By the time the effects of 
a reverse shock begin to appear in the line profiles (around 
day 5000), our models imply that the grains could already 
be as large as several microns in radius and are likely to 
be larger than ~ 0.6 /im. Grains as large as this are more 
likely to survive destruction by sputtering in supernova re¬ 
verse shocks and in interstellar shocks (Silvia et al. 2010, 
2012; Slavin et al. 2015). It has been suggested that very 
large grains (radii up to 4.2 /an) formed in the ejecta of SN 
2010jl within a few hundred days after the explosion (Gall 
et al. 2014). The grain radii that W15 and ourselves obtain 
for SN 1987A at very late epochs are nearly as large as found 
by Gall et al. (2014) for SN 2010jl, with both results sug¬ 
gesting that grains large enough to survive the destructive 


force of a reverse shock have formed by a few hundred days 
post-explosion. 

The dust masses obtained from our modelling of 
SN 1987A’s line profiles support the conclusion of W15 that 
even after ~3000 d the dust mass was still only a fraction of 
its current value. This contrasts with the results of Sarangi 
& Cherchneff (2015) whose grain chemistry models predict 
that ejecta dust masses should plateau by around 5 yr af¬ 
ter the explosion. Our results show that SN 1987A’s dust 
mass had reached of the order of 0.1A7q by day 3604. Since 
its present dust mass is several times larger than this (Mat- 
suura et al. 2015, W15), a substantial fraction of the current 
dust mass must have condensed after this epoch, in agree¬ 
ment with the conclusions of W15. 

Ideally, our models would cover the entire evolution of 
SN 1987A’s Ha line profiles up to the present day. How¬ 
ever, the excitation of gas in the outer edges of the ejecta 
by the reverse shock after ~ day 5000 results in significant 
broad and asymmetric emission that dominates the origi¬ 
nal line profile (Fransson et al. 2013). In addition to this, 
the narrow lines from the equatorial ring start to become so 
strong relative to the declining broad Ha profile that, post¬ 
removal, not enough of the broad profile remained to be 
able to reliably infer information from the profile structure. 
These factors may be common to some other CCSNe that 
have interactions with surrounding circumstellar material. 
Care should also be taken to ensure that any observed late¬ 
time line profiles being modelled are not in fact the product 
of a light echo reflecting the spectrum from near maximum 
light. Nonetheless, detailed line modelling of asymmetric line 
profiles has proved effective in determining dust masses in 
the ejecta of SN 1987A at multiple epochs during the first 
ten years after outburst. The method clearly has wider ap¬ 
plication to other supernovae. 


7 CONCLUSIONS 

We have investigated the effects of scattering and absorp¬ 
tion by ejecta dust on supernova line profile shapes and the 
different characteristic features that may be produced. In 
particular, attention is drawn to the fact that a classical 
blueshifted peak and asymmetric profile with most flux on 
the blue side is not the only profile type that can signify the 
presence of dust. In the case of strong dust scattering, line 
profiles can have the majority of their flux on the red side. 
Even with just some dust scattering, profiles can often ex¬ 
hibit an extended red scattering wing, although care should 
be taken to ascertain that this cannot be accounted for by 
electron scattering (electron scattering optical depths should 
usually only be significant at very early epochs, < 200 d). 
The line peak should always he on the blue side, with a line 
peak velocity that will often correspond to the minimum ve¬ 
locity at the inner edge of the ejecta shell. If not obscured by 
narrow circumstellar [N il] 6584 A emission, a pronounced 
shoulder or corner may be present on the red side of the 
profile, also corresponding to the minimum velocity at the 
inner edge of the ejecta shell. 

We have modelled the Ha and [O i] A6300,6363 A line 
profiles from SN 1987A over a range of epochs and have 
obtained dust masses of the order of 0.1 Mq by day 3604. 
We derive a sigmoid fit to our dust mass data that predicts 
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a current dust mass of 0.68 Mq, in line with current SED- 
based dust mass estimates for SN 1987A. We find that large 
grains are necessary in order to reproduce both the extended 
red scattering wings and the asymmetry seen in several of 
the lines and that grains larger than 0.6 /im have formed 
by day 714, while by day 3604 grain radii of ~ 3.5 /im are 
needed. We find from fits to the Ha profile that dust masses 
cannot have exceeded a fewxl0 _,i M@ on day 714 for all the 
grain types investigated, apart from glassy pure magnesium 
silicate grains, for which up to 0.07 M@ can be fitted. 

The observed red-blue line asymmetries persist right 
through to day 3604 and beyond - if no further dust had 
formed after day ~800 then the expansion of the ejecta 
dust shell would cause dust optical depths to drop rapidly 
with time thereafter, leading to the disappearance of red- 
blue asymmetries. Just to maintain the observed degree of 
red-blue asymmetry seen at the earlier epochs therefore re¬ 
quires that dust must have continued to form beyond those 
epochs. 
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that it is unnecessary to consider terms of the order of 0( i r) 
and thus A may be reduced to 
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The new direction of travel and frequency in the observer’s 
frame are therefore given by 

v = v(l - xp x - yP y - zp z ) (A4) 

x = Tp(x - p x ) 

V = ^(x~Py) 
z = Tp(x-fl z ) 

For each scattering event, the packet must be trans¬ 
formed both into and out of the comoving frame. The reverse 
transform is applied by using the inverse Lorentz matrix A -1 
which is obtained by reversing the sign of v. Positive v is 
defined for frames moving away from each other and thus v 
is defined to be negative in the direction of the observer. 


APPENDIX A: THE LORENTZ TRANSFORM 
FORMALISM 

Since the outflow velocities in supernovae are high, the pho¬ 
ton packets are subject to Doppler shifting upon emission 
and at each scattering event. When the packet is initially 
emitted, it has a frequency and a trajectory in the rest frame 
of the emitter. Both of these must be transformed to the 
observer’s frame in order for the packet to be propagated 
through the grid. The new direction and frequency in the 
observer’s frame may be simply found by transforming the 
momentum four-vector P which is defined as 


This paper has been typeset from a JpX/JATfrX file prepared by 
the author. 
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We may then derive P', the momentum 4-vector in the ob¬ 
server’s frame using the relation 


P' = A P (A2) 

where 
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In practice, the velocities considered are low enough 
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